Chinese Journal of Aeronautics, (2015), 28(5): 1343-1354

JOURNAL OF

AERONAUTICS

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

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

Autonomous gliding entry guidance with geographic c^mtk constraints

Guo Jie *, Wu Xuzhong, Tang Shengjing

Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China

Received 14 August 2014; revised 11 March 2015; accepted 25 May 2015 Available online 1 September 2015

KEYWORDS

Active disturbance rejection control;

Drag-versus-energy profile; Extended state observer; Geographic constraints; Guidance strategy; Trajectory generation

Abstract This paper presents a novel three-dimensional autonomous entry guidance for relatively high lift-to-drag ratio vehicles satisfying geographic constraints and other path constraints. The guidance is composed of onboard trajectory planning and robust trajectory tracking. For trajectory planning, a longitudinal sub-planner is introduced to generate a feasible drag-versus-energy profile by using the interpolation between upper boundary and lower boundary of entry corridor to get the desired trajectory length. The associated magnitude of the bank angle can be specified by drag profile, while the sign of bank angle is determined by lateral sub-planner. Two-reverse mode is utilized to satisfy waypoint constraints and dynamic heading error corridor is utilized to satisfy no-fly zone constraints. The longitudinal and lateral sub-planners are iteratively employed until all of the path constraints are satisfied. For trajectory tracking, a novel tracking law based on the active disturbance rejection control is introduced. Finally, adaptability tests and Monte Carlo simulations of the entry guidance approach are performed. Results show that the proposed entry guidance approach can adapt to different entry missions and is able to make the vehicle reach the prescribed target point precisely in spite of geographic constraints.

© 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

As a primary focus of entry technology, entry guidance is concerned with steering the vehicle to the designated target with prescribed conditions while satisfying all necessary path

* Corresponding author. Tel.: +86 10 68918678. E-mail address: southcross@yeah.net (J. Guo). Peer review under responsibility of Editorial Committee of CJA.

constraints in spite of the internal and external disturbances of system. Path constraints are usually designed for the sake of vehicle's safety and operation and the typical ones include heating rate, aerodynamic load and dynamic pressure which would present a great challenge for entry guidance design. Moreover, when it comes to the relatively high lift-to-drag ratio L/D entry vehicles, the large longitudinal phugoid oscillations of the gliding trajectory further increase the difficulty in designing effective flight control systems which are usually based on aerodynamic control effectors. As a result, large longitudinal phugoid oscillations are needed eliminating.1 Furthermore, the relatively high L/D vehicles possess a capability of lateral maneuvers without power. Additional flight

http://dx.doi.org/10.1016/j.cja.2015.07.009

1000-9361 © 2015 The Authors. Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

mission requirements such as flying over specified waypoints and avoiding prescribed no-fly zones may be introduced to the long-range maneuverable hypersonic entry vehicles. To be more specific, waypoints are predetermined positions for multiple payload deployments or reconnaissance missions, and no-fly zones are exclusion zones that cannot be passed for threat avoidance or due to geopolitical problem.2 The way-points and no-fly zones are collectively referred to geographic constraints. Jorris and Cobb3'4 first optimized the entry trajectory that satisfied geographic constraints for the hypersonic cruise vehicle (HCV) in a two dimensional platform and extended the approach into the three-dimensional models for common aero vehicle (CAV). Zhao and Zhou5 presented a multi-phase technique based on Gauss pseudo-spectral method to generate an optimal reentry trajectory satisfying geographic constraints. The development of pseudo-spectral method makes it possible to optimize the entry trajectory subject to typical path constraints onboard.6-8 Furthermore, some novel fast trajectory planning methods are developed from the previous idea.8'9 Xie et al.2'10 proposed two trajectory planning methods considering geographic constraints based on the equilibrium glide assumption and drag-versus-energy profile, respectively. However, both of the planning methods have not been extended to entry guidance algorithm.

Generally, the entry guidance algorithm can be grouped into two categories: guidance using predictor-corrector and guidance using reference trajectory.11 The predictor-corrector guidance predicts the terminal condition employing analytical or numerical propagations and adjusts the design parameters so that the errors of terminal constraints can be corrected. The guidance algorithm enjoys the main advantage of little need for a pre-stored reference trajectory but enlarge the difficulty for predictor-corrector guidance in enforcing various path constraints strictly.12,13 Joshi et al.14 proposed a predictor-corrector guidance considering typical path constraints. Moreover, for the aforementioned high L=D entry vehicles, a unique phenomenon is that the entry trajectory is prone to large oscillations. To eliminate those oscillations, Xu et al.15 combined a predictor-corrector guidance with Quasi-Equilibrium glide condition (QEGC) which made the generated trajectory smooth, but also enabled the guidance to enforce typical path constraints. Lu et al.1,16 presented an effective feedback compensation in conjunction with the predictor-corrector guidance algorithm to eliminate phugoid oscillations.

For guidance using reference trajectory, guidance algorithm comprises a trajectory planner and a trajectory tracker. Shuttle entry guidance17,18 uses drag acceleration as a surrogate control variable. The drag acceleration profile can be mostly designed offline according to trajectory optimization algorithm, while minor online adjustments are also utilized to mollify predicted down-range errors. The required downrange is obtained under the assumption that the vehicle flights on a great circle arc to the target point. The drag tracking law specifies the magnitude of a bank angle, while the sign of bank angle is determined by a heading error corridor. With great success in application, shuttle entry guidance has become a baseline approach for many entry vehicles. Mease et al.19 presented a reduced-order entry trajectory planning method which was also based on drag acceleration profile. Since the strategy combines the longitudinal and lateral motion, the trajectory planning methods realize large crossrange entries. Saraf

et al.20 developed Mease's method and presented the design and performance assessment of evolved acceleration guidance logic for entry (EAGLE). Leavitt and Mease21 presented trajectory planner by constructing a more easily tracked drag profile using interpolation, meanwhile the planner possessed near-maximum downrange and crossrange capabilities. With the aid of the QEGC, Shen and Lu22,23 introduced the altitude-versus-velocity profile and divided entry trajectory into initial descent phase, quasi-equilibrium glide phase and pre-TAEM phase. Li et al.24 proposed a unified analytical expression to describe the three typical path constraints, which greatly simplified the reference trajectory design.

The trajectory tracker follows the reference trajectory profile during flight. Shuttle entry guidance employs a gain-scheduled proportional-integral-derivative (PID) control logic to track the reference trajectory.17 However, PID is problematic for tedious gain scheduling and tuning issues. As a result, Mease25 and Lu26 et al. introduced guidance laws based on feedback linearization and Morio et al.27 proposed general guidelines for the application of flatness theory to the entry guidance. The trajectory tracker can also be designed by optimal control theory. Dukeman28 utilized linear quadratic regulator theory to track the reference profile. This method treated the trajectory-tracking problem as a regulation problem in the state space about the reference trajectory. Since the online computational burden associated with solving Riccati differential equation is an impediment for applications of the method, Lu29 introduced a method that is based on the approximate state and control at discrete time points to avoid the online integration of the Riccati differential equation. Similarly, the indirect Legendre pseudo-spectral method30,31 and the generating function method32,33 are also efficient numerical algorithms capable of freeing the equation solving from online integration.

Active disturbance rejection control (ADRC) proposed by Han34 is an emerging nonlinear control algorithm. ADRC inherits the quality of PID but combines with nonlinear control strategy. Vincent's tests35 show that ADRC control makes a significant improvement over the PID control under the same conditions. ADRC control method employs an extended state observer (ESO) that can precisely estimate the internal and external disturbances of the system and dynamically compensate the system accordingly. That makes the control method not depend on the accurate mathematical of the unknown object model. Compared to PID, ADRC outstands with more static and dynamic performance, strong robustness and adaptability.

The overall objective of this study is to develop a novel three-dimensional autonomous entry guidance considering the geographic constraints. In this paper, the feasible reference trajectory is designed onboard based on drag-versus-energy plane and the trajectory tracking method is designed based on ADRC. The other parts of this paper are organized as follows. Section 2 describes the entry dynamics and the constraints for the entry guidance problem. Section 3 presents the strategy of the rapid entry trajectory planner. Then a novel tracking law for altitude tracking based on ADRC is presented in Section 4. In Section 5, the performance of the proposed guidance method is assessed by numerical adaptability tests and Monte Carlo simulations. Finally, the conclusions are presented in Section 6.

2. Preliminary

2.1. Translational equations of motion

A group of nonlinear differential equations governs entry dynamics. With the time t as an independent variable, the dimensionless form of the three-dimensional (3-D) point mass dynamic equation of an entry vehicle over a spherical rotating Earth can be written as

■ dx ft \ x = dt = f(X; u)

where the state vector x — [r, 0, /, V, y, W] and the control vector u — [a, r]. r is the radial distance from the Earth center to the vehicle, 0 the longitude, / the latitude and V the magnitude of Earth relative velocity. The flight path angle y is the angle between the Earth relative velocity vector and the local horizontal plane, and the heading angle W is the angle between the projection of same velocity vector on the local horizon plan and local longitude line, measured from the north in the clockwise direction. Moreover, the control variables used are angle of attack a and bank angle r, the latter is defined as the angle from the vertical upward direction to the lift vector, positive to right.

In dimensionless form, time, length, velocity, acceleration, angular rate, mass and density are normalized by tS, lS, VS, aS, wS, mS and qS, respectively. These dimensionless quantity is modeled using Eq. (2):

fts, Is, Vs, aS, «s, mS, Ps}

— { \/re/go, re, Pgr", go, VgO/Te, mv

where mV is the mass of vehicle, re the equatorial radius of Earth, g0 — i/r2 and i the planet's gravitational constant. The non-conventional energy E is defined as10

E — V2/2 - 1 /r (3)

If the Earth-rotation term of V in Eq. (1) is ignored, we can

, mv/r3}

dE/dt = -VD

where D describes the dimensionless form of the drag acceleration. The non-conventional energy E decreases along with the trajectory given the assumption D p 0 and V P 0 in the definition. Since time is not typically a critical parameter in guidance problem, the monotonically decreasing variable E should be treated as an appropriate independent variable of the entry dynamic equations for trajectory planning.

The dimensionless form of 3-D equations of Eq. (1) can be rewritten with E as independent variable as follows20: sin y

/' = -

Dr cos / cos c cos w

, i l • COer . , . ,

V = — — 2 sin c cos /(sin c cos / — cos c cos w sin <

V VDr2

■ cos r +

V2D C e2 r cos / V2D cos c

cos c 2C e

VD — VD

sin w cos /

(cos c cos / + sin c cos w sin (

W' = —

L cos c sin W tan / 2œe .

- sin r----(sin/

V2D cos c Dr VD

— tan c cos W cos /) — 2

V2D cos c

sin w sin / cos /

where (■)' denotes d(-)/dE, Me describes the dimensionless form of the angular rate of the Earth-rotation, L means the dimensionless form of the lift acceleration, and lift is rotated by the bank angle. Note that V can also be obtained by Eq. (3) without integrating over non-conventional energy of Eq. (8). The control variable r appears explicitly in the equations of motion, whereas another control variable a appears through D and L. The drag acceleration D and the lift acceleration L take the form of

D = 1 qV2 — Cd (a, Ma)

L = 1 qV2 Sref CL(a, Ma) 2m

where Sref is the dimensionless form of the reference area and m refers to the dimensionless form of the mass of the vehicle, CD and CL stand for drag and lift coefficients, respectively. Both of the coefficients depend on the angle of attack a and (2) Mach number Ma. In this paper, CD and CL are obtained by36

CD = 0.16721 — 0.01853Ma + 0.00151a

+0.00074Ma2 — 0.00045Maa + 0.00081a2

CL = —0.04702 — 0.00673Ma + 0.05071a

+0.00042Ma2 — 0.00096Maa + 0.00043a2

q is the dimensionless form of the atmospheric density that can be given according to an exponential function of altitude

q = q0 exP

where q0 refers to the dimensionless form of density at the reference radius re, and hs is the dimensionless form of the constant scale height.

2.2. Crossrange and downrange

As shown in Fig. 1, O(00, /0) is the intersection of the Earth's surface and the connecting line that connects entry interface (EI) and the center of the Earth. T(0f, /f) is the target point and E(0e, /e) is the current position of the vehicle. o0f is the angle between the tangent of great circle arc OT and the line of longitude at O. o0e is the angle between the tangent of great circle arc OE and the line of longitude at O. Both o0f and o0e are measured from north in clockwise direction. The downrange RD and the crossrange RC for vehicle's current state can be obtained by

sin Rc — sin ^0e sin(t0f - %;) (16)

cos Rd — cos b0e/cos RC (17)

where b0e is the length of the great circle arc OE.

Fig. 1 Schematic map of downrange and crossrange. 2.3. Path and terminal constraints

The upper limits of the stagnation point heating rate QQ, the dynamic pressure q and the aerodynamic load n are given by

Q = cqqp0-5V3-l5(mV5g0'575r0-075) 6 QQm q = 0 . 5q F2 (mvg0r;2) 6 qmax n =\JL2 + D2 6 nmax

where cq is a vehicle dependent constant, QQmax, qmax, and nmax are the maximum accept values of heating rate, dynamic pressure and aerodynamic load, respectively. All of the three typical path constraints above are considered as rigid constraints.

The flight path angle y is small and varies slowly for the major portion of gliding entry trajectory. Assuming cos y = 1, y' = 0 and r = 0 in Eq. (9) and ignoring the Earth-rotation term, we can obtain

1 V2 L - - + — P 0

The entry vehicle does not pose a risk if the inequality constraint is invalid, which serves to reduce the phugoid oscillations in altitudes along with the entry trajectory.22 Hence, the constraint does not need enforcing strictly and the zero-bank QEGC is usually employed as a soft constraint. Note that the assumptions cos y — 1 and y' — 0 are only employed to deduce QEGC.

The geographic constraints comprise waypoints constraints and no-fly zone constraints. Waypoints are the predetermined positions for the missions of multiple payload delivery or reconnaissance. In this study, the longitude and latitude of the waypoints are constraints; however, time, control variables, altitude, velocity, flight path angle and heading angle are not considered as constraints. If B(hb, /b) is the point at entry trajectory that is at the minimum distance from the way-point A(ha,/a), the constraint for the waypoint A can be expressed as

|hb - ha| 6 8g, |/b - /aI 6 8/ (22)

where eh and e/ are the preselected small positive values.

No-fly zones are the exclusion zones that cannot be passed for threat avoidance or due to geopolitical restrictions. In this study, the no-fly zones are specified as the circular exclusion zones with infinite altitudes. If the central position of one no-fly zone is (hc, /c) and the radius is rc, the constraint for the no-fly zone can be expressed as

(h - hc)2 + (/ -

/c)2 > r2

where (h, /) is a point at the entry trajectory. Note that the unit of rc is in radians.

The typical final condition for an entry trajectory is that the trajectory reaches a position at a specified distance sf .1 In other words, the terminal radial distance r(tf), the terminal velocity V(tf) and the great-circle range to the target point s(tf) are specified by

r(tf) = r? V(tf) = Vf

s(tf)=s(0f', /f) = sf

where the great-circle range s refers to a function of longitude and latitude.

Control variables and angular rates of control variables are considered as constraint in this study. a, r, a and r are all restricted within a certain range of values. Therefore, the control variables can be modeled as

a e [amin, amax], r e [rmin, ГшиРI 6 IaImax, I<TI 6 I<TIm

where the subscripts ''min" and ''max" denote the minimum and maximum acceptable values, respectively.

3. Trajectory planning approach

The boundaries of entry corridor can be determined uniquely once the angle of attack a is specified. In this study, the nominal angle of attack is scheduled as piecewise linear functions of velocity and can be expressed as

■■={ al + V(V- Vl)

V P V1

V1 > V > V2

where Fb V2, ai, are the parameters of the angle of attack profile. The vehicle dependent entry corridor in drag-versus-energy plane is designed offline. With the determination of nominal angle of attack by Eq. (28), the inequality constraints of Eqs. (18)-(20) can be converted into drag acceleration constraints. The maximum and minimum boundaries of the D-E corridor Dmax(E) and Dmin(E) are given by

Dmax(E) = minfDg (E), Dq(E), Dn(E)}

Dmin(E) = Dqegc(E)

where DQ_ , Dq

and Dn represent the drag acceleration constraints corresponding to the maximum heating rate, dynamic pressure and aerodynamic load, respectively. Dqegc describes the drag acceleration constraints for zero-bank QEGC.

Since this study does not address the problem of selecting a feasible entry target, it assumes that a solution to the trajectory-planning problem exists. The objective of the

proposed trajectory planner is to generate a feasible trajectory onboard that meets the boundary conditions and the path constraints and also fulfills the equations of motion with the nominal models.

The entry trajectory is divided into two phases: the initial descent phase and the glide phase. The trajectory in each phase features distinctive characteristics. For the initial descent phase, QEGC is not valid and the vehicle rarely possesses lateral mobility. For glide phase, the atmospheric density becomes relatively high and the lift is sufficient for lateral maneuver that could satisfy the geographic constraints. The interface between the two phases is named as the transition point.

3.1. Initial descent phase

In the initial phase of entry flight, the zero bank QEGC is not valid because the atmospheric density is too low such that the aerodynamic lift cannot afford the weight of vehicle to maintain its flight path angle. Starting from the EI, the glide entry vehicle needs to descend and enter the entry corridor. The nominal angle of attack profile modeled using Eq. (28) and a suitable constant bank angle r0 are used to integrate non-conventional energy of Eqs. (5)-(10). The appropriate magnitude of r0 is determined by starting from zero and increasing the magnitude by a fixed increment progressively until both the flight constraints and the transition conditions are satisfied. The transition point between the initial and glide phases takes the form

|dr/dV - (dr/d V)qegc I 6 «QEGC

where £qegc is a preselect small positive value, dr/dV — — Vsiny(D + gsiny), and (dr/dV)qegc is the slope of the QEGC and the current point (r, V) and it can be obtained by differentiating the QEGC with respect to V.

The sign of initial bank angle r0 is determined as follows. As shown in Fig. 2(a), if the first geographic constraint is a waypoint, the initial sign of bank angle can be specified as

sign(ff0) = -sign(W0 - twP1 )

where W0 is the initial heading angle and oWP1 is the angle between the connecting line and the line of longitude at EI. The connecting line refers to the line that connects EI and the first waypoint. As shown in Fig. 2(b), if the first geographic constraint refers to the no-fly zone, the initial sign of bank angle can be modeled as

sign(ff0) = -sign(W0 - tA1 )

where oA1 refers to the angle determined by the tangent line of the first no-fly circle zone intersecting the longitude line at EI. In addition, oA1 is measured from the north in clockwise direction.

3.2. Glide phase

To facilitate the presentation, it assumes that the energy is normalized as

E - E0

Ef — E0

Fig. 2 Relative position between EI and first geographic constraint.

In the above equation, E is the non-conventional energy of the current state. E0 and Ef stand for the initial non-conventional energy and terminal non-conventional energy, respectively. Thus E equals 0 in the initial state and equals 1 in the state of the target point.

3.2.1. Longitudinal sub-planner

A longitudinal sub-planner is designed to specify the drag-versus-energy profile online to meet the path and terminal constraints including trajectory length, terminal altitude, terminal magnitude of velocity, heating rate, aerodynamic load, dynamic pressure and QEGC. One of the most significant objectives in the drag planning is to achieve the downrange of a target point. If the reference drag acceleration profile Dref is specified, the trajectory length S can be obtained by

r 1 -dE (34)

Dref(E)

Motivated by the computing time, we employ interpolation to obtain the reference drag profile. Ref. 21 indicates that the trajectory generated from drag interpolation is feasible to sufficient accuracy if both the lower and upper drag profiles are associated with feasible trajectories. A weighted combination of minimum and maximum drag profiles is used to generate the drag acceleration profile. In this study, the interpolation depends on the parameter k. If k — 1, the minimum drag profile is generated, and if k — 0, the maximum drag profile is generated. The reference drag profile Dref is given by

Dref(E) = kD max (E) + (1 - k)Dn

n(E) (35)

in are obtained by Eq. (29)

where 0 < k < 1, Dmax and . offline.

Since the drag acceleration of both the transition point and the target point is specified for the trajectory planning of glide phase, we use two linear functions with the energy as abscissa to modify the reference drag profile:

Di - Dc Ei — Ec

Dref = <

kDmaK +(1 - k)Dm

D2 - Df E2 — Ef

(E - Ec) Ec 6 E < Ec + 0.1 Ec + 0.1 6 E < 0.9

(E - Ef) 0.9 6 E < 1

where Dc and Ec are the drag acceleration and the non-conventional energy of the transition point, while Df and Ef are those of the target point. D1 and E1 are the drag acceleration and the non-conventional energy for E — Ec + 0.1, D2 and E2 are the drag acceleration and the non-conventional energy for E — 0.9. After Eqs. (34) and (36) are combined, the trajectory length only depends on the parameter k.

The reference magnitude of the velocity Vref and the reference radius rref (the reference altitude href) can be obtained by solving the set of Eqs. (3) and (11) in the case that the reference drag-versus-energy profile is specified. The equation set is solved by secant method in this paper. Then the reference flight path angle yref can be obtained by

Cref = arcsin( r'refDref)

where /ref can be obtained by finite difference method. Ignoring the Earth-rotation terms of Eq. (9), the reference bank angle rref can be expressed as

VefAef ( , , ( l V2\ Cosf

cos rref = -JT- ( "Cref + — -— ) T7TT, ) (38)

ViefDref.

where yref can also be obtained by finite difference method and Lref is the lift acceleration corresponding to the Dref.

3.2.2. Lateral sub-planner

The constraints of the waypoints and the no-fly zones are considered in lateral trajectory planning, and satisfied by reversals of bank angle. The reduced-order equations of lateral motion for planning frame derives from the equations of motion in 3D Eqs. (5)-(10) with the non-conventional energy E as an independent variable.

w = ■

VrefDref cos Cref

cos yref sin W tan /

Drefrref

VrefDr

- (sin / - tan yref cos W cos <

Wi^ref

/ = ■

VrefDref cos Cref cos Crefsin W

Dfref cos / cos Cref cos W

Drefrref

sin W sin / cos /

where Lref, Dref, rref, Vref, cref and magnitude of the bank angle r can be determined by longitudinal sub-planner. Hence Eqs. (39)-(41) can be solved if the sign of bank angle is specified.

(1) Reversal logic-based on two-reverse mode for waypoints and target point.

A series of points including the transition point, waypoints and the target point, are referred to as the reference points (RP). The longitude and latitude information of RP can be transformed into the downrange and crossrange information

by Eqs. (16) and (17). Hence the ,th RP can be expressed as RP,-(Rd,-, Rc,) instead of RP,(0rp,-, /rk), where RDi, Rqu Örp,-, /RPi are downrange, crossrange, longitude and latitude of RP,-, respectively. The crossrange errors are expected to be zero at both the waypoints and the target point. Thus the trajectory plan approach for them can be designed in a similar way. In the process of entry, we can predict the crossrange at the way-points and the target point by integrating the equations of lateral motion Eqs. (39)-(41).

If one reverse mode is used between two adjacent points RP,_1 and RP,-, the parameter k cannot be updated after the reverse point and the reference trajectory cannot, either. Hence a single bank reversal can hardly keep a prefect control effect, and the crossrange error at RPi may be intolerable if the reverse point is too far away from RP,-. On the other hand, unnecessary reversals are not expected either. Hence we use two-reverse mode for the crossrange error requirement at RP,. We define the dRD as

dRn, = Rd

Therefore, the down ranges of two reverse points are specified as Rd,-1 + P1dRDi and Rd,-1 + P2dRDi, where 0 < p1 < p2 < 1 are the parameters of the two reversal mode. The reversal strategy is as follows:

(A) If the downrange of current state fulfills RDi_1 6 Rd < RDi_1 + p1dRDi, p2 is specified as a constant value close to 1.0, and the lateral sub-planner updates p1 to minimize the crossrange error at RPi.

(B) If the downrange of current state RDi_1 + p1dRDi 6 Rd < RDi_1 + p2dRDi, p2 is not specified any more, and the lateral sub-planner updates p2 to minimize the crossrange error at RP,.

(C) If the downrange of current state can be described as RDi-1 + p2dRDi 6 Rd < RDi, the lateral sub-planner does not update either p1 or p2.

(2) Reversal logic based on dynamic heading error corridor for no-fly zone avoidance.

The reversal logic based on dynamic heading error corridor is an extension of shuttle-type lateral guidance logic. It is essential for defining dynamic heading error corridor to determine the minimum and maximum boundaries of the head angle W according to the relative motion between no-fly zone and vehicle current line of sight. As shown in Fig. 3, we assume that the central position and the radius of the no-fly circular zone are N(hn, /n) and rn, and the vehicle's current and target position are £(he, /e) and T(hf, /f). tef is defined as the angle between the vehicle's current line of sight ET and local longitude line. oeg is defined as the angle between the tangent EG and local longitude line. Both oef and oeg are measured from the north in the clockwise direction.

As shown in Fig. 3, the vehicle is on the west of the target point and avoids the no-fly zone from north side. There are two possible cases for relative position of line of sight ET and no-fly zone. One is that ET overpasses the no-fly zone (Fig. 3(a)) and the other is that ET does not overpasses the no-fly zone (Fig. 3(b)). If ET overpasses the no-fly zone, the boundaries of the heading angle are specified as

(b) ETdose not overpass no-fly zone Fig. 3 Position of line of sight relative to no-fly zone.

Wmin — °eg — Wmax — °eg

If ET does not overpasses the no-fly zone, the boundaries of the heading angle are given by

Wmin — min(teg — Oef — d//2) Wmax — min(teg, ^ef + d//2)

According to the previous analysis, the overall approach to determine the range of the heading angle is as follows. The parameters s1, s2 and s0 are defined as

1 he P hf

— 1 he < hf

1 /e P /n (45)

— 1 /e < /n

s0 — s:s2

where s1 — 1 means that the vehicle is on the east of the current target point and si — — 1 means the vehicle is on the west of the target point. In addition, s2 — 1 happens when the vehicle avoids the no-fly zone from the north side, while s2 — 1 happens when the vehicle avoids the no-fly zone from the south side.

(A) If s1 he > s1hn — rn/re, the no-fly zone has not been avoided. The minimum and maximum values of the heading angle are given by

i Wmin — s0 max (s0 min(teg, teg + s0dW), s0 (tef — ¿W/2)) \Wmax — s0 max(s0max(teg, teg + s0dW), s0 (tef + dW/2))

(B) If s1 he 6 s1hn — rn/re, the no-fly zone has been avoided and the values of Wmin and Wmax are

Wmin — tef - d//2 Wmax — °ef + d//2

Since the minimum and maximum boundaries of the head angle W are specified, the reverse logic based on dynamic heading error corridor is designed as

{-1 W > Wmax

1 W < Wmin (48)

sign(r,-1) Wmin 6 W 6 Wmax

where C7_i is the latest bank angle command.

In summary, the main function of the longitudinal sub-planner is to adjust the parameter k to meet the terminal down-range requirement and the lateral sub-planner specifies the appropriate bank-reversal timings so that the vehicle can fly over the waypoints, avoid the no-fly zones and meet the terminal crossrange requirement. The longitudinal and lateral sub-planners are iteratively employed until the downrange requirements, crossrange requirements and all the path constraints are satisfied. The flowchart of the proposed trajectory planning approach is shown in Fig. 4. The iterations are not more than 3 in the following test in Section 5, hence the trajectory planner can generate the feasible glide phase trajectory rapidly.

4. Trajectory tracker based on ADRC

According to Section 3.2, radius r, velocity V, heading angle W and coordinates (h, /) of the transition point are specified for the trajectory planning algorithm at glide phase. However, the flight path angle of the transition point is not specified, which makes the vehicle deviate from the reference trajectory when the reference bank angle rref is employed. Moreover, the approximation in Eq. (4) also results in deviation. To solve the problem, a trajectory tracker based on ADRC is introduced in this section. The objective of the trajectory tracker is to steer the bank angle so that the vehicle follows the reference trajectory produced by the planner. As shown in Fig. 3, the controller, which has the function of ADRC, consists of following parts34: (A) a nonlinear tracking differentiator (TD) that is used to arrange the ideal transient process of the system; (B) ESO, which estimates all the disturbances from the system output and compensates the disturbances according to the estimated values; (C) a nonlinear state error feedback (NLSEF) that is used to obtain the control of the system.

In Fig. 5, V stands for reference signal. Altitude h is selected to be reference signal in this study. The altitude rate with degree g — 2 can be obtained as

h = a + bu

Note that u — cos r is the control variable. In addition, a and b can be expressed as

V2 i 2 a — — D siny H--cos y —j + «2r sin y cos/(sin y cos /

— cos y cos W sin /) + 2we V cos y sin W cos / + wjr x cos/(cos y cos / + sinycos W sin/) (50)

b — L cos y

and the altitude rate can be obtained as h — V sin y

( End J Fig. 4 Flowchart of trajectory planning approach.

Define the state x1 — h, x2 — h, and x3 — a, where, x3 is the extended state. According to Eq. (49), the system can be rewritten as

x 1 — x2 x2 — x3 + bu x 3 — v

.y — x1

where v stands for the derivative of a. Then, an extended state observer (ESO) for Eq. (53) can be obtained as

e — Z1 — y Z — Z2 — fte z2 — z3 — b2fal(e, d, d) -z3 — —b3fal(e, d1, d)

where e is the estimation error of ESO, and z1, z2, and z3 are the observer outputs. z1 is used to estimate system output, z2 is used to estimate the differential of system output and z3 is extended state variable to estimate comprehensive disturbance. b1, b2, and b3 are the observer gains. fal(-) is a continuous power function with a linear segment near origin and can be expressed as

fal(e, d, d) —

|e| 6 d > d

where 0 < d < 1 , d > 0.

d, d1, d, b1, b2, and b3 are the parameters to be designed. The output of the observer can track h, h, and a precisely if the six parameters are tuned appropriately. Using the ADRC algorithm, the control law can be designed as

u —(% — z3)/b0 (56)

where b0 is the approximate value ofb. u0 can be obtained by a nonlinear combination of error signal and its differential:

u0 — b01fal(e1, a1, d) +b02fal(e2, a2, d) (57)

where the conditions are 0 < a1 < 1, 0 < a2 < 1 and e1 — h — z1, e2 — h — z2.

Furthermore, the angular rate of the bank angle especially at the reverse point is also considered and the bank angle can be modified by

( sign(rref)|c

r — < sign(rref)|r| — r— 1

r;—1 + ^-;-—;-T r m

|sign(rref )|r| — ri—11

|r | 6 rm |<r | > rm

where rref is the reference bank angle generated by trajectory planner. r and rl—1 are current and latest bank angle command, respectively. r and rl—1 are generated by trajectory tracker.

5. Assessment of guidance performances

As stated above, the guidance is composed of a trajectory planner and a trajectory tracker. The proposed trajectory planner updates the reference trajectory every 60 0 s, however, as the vehicle flies close to the target, that is, the normalized energy E P 0 . 8, the update period should be reduced and the author specified it as 20 0 s. The reference trajectory generated by trajectory planning each time is employed as the initial guess for next time. In this paper, the initial guess of parameter k in Eq. (36) is 0.5. In addition, the initial guess of parameters p1 and p2 are 0.5 and 0.95, respectively. The preselected values

Fig. 5 ADRC topology.

Table 1 Positions of EI and geographic constraints for 5 missions.

Mission Initial longitude and latitude (h '0> /0)(°) Initial heading angle W0(°) Geographic constraint

Mission 1 (55.0, -20.0) 65.0 None

Mission 2 (55.0, -10.0) 65.0 Waypoint A

Mission 3 (55.0, 0) 80.0 Waypoint B, Waypoint A

Mission 4 (55.0, 10.0) 90.0 No-fly zone C, Waypoint A

Mission 5 (55.0, 20.0) 85.0 No-fly zone C

of e0 and e/ in Eq. (22) are set as 0.15 and 0.2, respectively. The independent variable of the reference trajectory is non-conventional energy E, while the practical trajectory is obtained by integrating over the time of the Eq. (1). After the reference trajectory is generated, the trajectory tracker based on ADRC generates the control command of the practical trajectory every 0.5 s.

The Lockheed-Martin CAV-H model is employed in the numerical tests. The weight of CAV-H model is about 907.2 kg and the lift and drag coefficients are obtained from Ref. 36 with the reference area prescribed as 0.4839 m2. The maximum lift-to-drag ratio is 3.5. The control constraints are a 2 [5°, 25°], r 2 [-90°, 90°], |a| 6 5(°)/s and |r| 6 20(°)/s. The parameters of the angle of attack profile Eq. (28) are specified as V! — 5000 m/s, V2 — 2000 m/s, a! — 20°, and a2 — 10°. In addition, the maximum heating rate, aerodynamic load and dynamic pressure are specified as QQmax — 1.5 x 106 W/m2, nmax — 3.0, qmax — 105 Pa.

5.1. Guidance adaptability assessment with different EI and different geographic constraints

The altitude, velocity and flight path angle of EI are specified as h0 — 80 km, V0 — 7000 m/s and y0 — 0°. The terminal constraints are h, — 27 km, V, — 2000 m/s, 0,5 — 115° and

— 15°.

The predetermined geographic constraints include waypoint A, waypoint B and no-fly zone C. The longitude and the latitude of waypoint A and waypoint B are (100.0°, 10.0°) and (85.0°, 5.0°), respectively. The central position of the no-fly zone C is (80.0°, 15.0°) and the radius of the no-fly zone C is 667.9 km.

In order to test the adaptability of the guidance method, 5 missions with various EI and geographic constraints are employed. The initial longitude, initial latitude, initial heading

201_,_,_,_._,_,_._,_

0 600 1200 1800 Time (s)

Fig. 6 Altitude histories of entry vehicle for Mission 1-Mission 5.

angle and geographic constraints are shown in Table 1. In addition, the results for the 5 missions without disturbances are shown in Figs. 6-8.

The altitude history shown in Fig. 6 illustrates that the large phugoid oscillations are effectively damped out in the entry trajectory. All entry trajectories are divided into initial descent phases and glide phases. For initial descent phase, the vehicle descends and enters the entry corridor. For glide phase, the reference altitude which is specified by the major part of the reference drag profile decreases monotonically, and the altitude tracker based on ADRC works effectively. As shown in Fig. 7, all the geographic constraints specified in Table 1 are satisfied and all the generated trajectories reach the same target point lastly.

The history of heating rate, dynamic pressure, aerodynamic load and drag acceleration is shown in Fig. 8. The bold lines in bottom right of Fig. 8 stand for maximum drag boundary and

Fig. 7 Ground tracks of entry vehicle for Mission 1-Mission 5.

Fig. 9 Altitude and flight path angle.

minimum drag boundary. It is observed that all the typical path constraints are satisfied during entry flight.

Furthermore, the test results compared with PID method for Mission 4 are shown in Figs. 9 and 10 based on the tracking law proposed by Section 4. In order to focus on the comparison of ADRC and PID, the reference trajectory, which trajectory planner generates for the first time, does not update.

As shown in Fig. 9, the ADRC altitude profile tracks the reference profile perfectly. In addition, the terminal altitude and velocity by the tracker based on ADRC are 27 009 km and 2000 .275 m/s, while PID's are 27 .058 km and 2002. 380 m/s. The flight path angle profile illustrates that the ripple and overshoot are greatly reduced with the help of ADRC controller. As shown in Fig. 10, the bank angle commands generated by PID and ADRC all satisfy the control constraints.

Fig. 8 Path constraints.

Normalized energy Fig. 10 Bank angle.

Fig. 12 Terminal altitude errors and velocity errors.

Table 2 Dispersions in Monte Carlo simulation.

Parameter Mean 3r

Initial altitude error (km) 0 1.0

Initial longitude error (°) 0 0.1

Initial latitude error (°) 0 0.1

Initial velocity error (km/s) 0 0.1

Initial flight path angle error (°) 0 0.1

Initial heading angle error (°) 0 0.1

Drag coefficient error (%) 0 10

Lift coefficient error (%) 0 10

Atmospheric density error (%) 0 10

Vehicle mass error (%) 0 5

60 70 80 90 100 110 120 Longitude (°)

Fig. 13 Ground tracks for 1000 dispersed cases.

Normalized energy n4 95 n4 9? n4 99 n5 01 1]5 03 ]15 05

Fig. 11 Altitude histories for 1000 dispersed cases. Longitude (°)

Fig. 14 Terminal positions for 1000 dispersed cases.

5.2. Monte Carlo simulation

0.130 km and 0.006 km, respectively. And the maximum and To further evaluate the robustness of the entry guidance mean values of terminal velocity error are 2.431 m/s and method proposed in this paper, dispersion study was executed 1.239 m/s, respectively.

for Mission 4 by Monte Carlo simulation. The distribution of As shown in Fig. 13, the ground tracks of Mission 4 for

the initial condition and modeling errors are shown in Table 2. 1000 dispersed cases avoid the no-fly zone C, then pass the The altitude history for 1000 dispersed cases is shown in waypoint A and reach target point successfully in the end. Fig. 11 and the zoom-in view near the target is shown in the The two zoom-in views near the waypoint A and target point inset at the top right of the figure. As shown in Fig. 12, are shown in the inset at the top and the bottom right of the the maximum and mean value of terminal altitude error are figure, respectively. Terminal positions for 1000 dispersed

cases are shown in Fig. 14. The maximum and mean values of terminal distance error are 5 .170 km and 2.902 km, respectively. It can be observed that 99.0% of the terminal miss distances are within 5 0 km of the target.

6. Conclusion

A framework of entry guidance methodology for relatively high L/D vehicles combining onboard feasible trajectory planner and robust trajectory tracker is established. The trajectory planner that is designed by drag-versus-energy profile can plan a feasible reference trajectory satisfying the geographic constraints online. Furthermore, the trajectory planning method is able to generate the associated reference altitude profile simultaneously. Large phugoid oscillations of the altitude profile are effectively damped out. Moreover, the robust trajectory tracker based on ADRC, is employed to track the reference altitude profile. The numerical adaptability test illustrates that the proposed guidance method features the capability of steering the entry vehicle to the target point precisely with different geographic constraints from different EI. Furthermore, the Monte Carlo simulation tests indicate that the proposed guidance method is robust to the presence of a wide range of initial uncertainties and dynamic disturbances.

Acknowledgement

This study was supported by National Natural Science Foundation of China (No. 11202024).

References

1. Lu P. Entry guidance: a unified method. J Guidance Control Dyn 2014;37(3):713-28.

2. Xie Y, Liu LH, Liu J, Tang GJ, Zheng W. Rapid generation of entry trajectories with waypoint and no-fly zone constraints. Acta Astronaut 2012;77(8-9):167-81.

3. Jorris TR, Cobb RG. Multiple method 2-D trajectory optimization satisfying waypoints and no-fly zone constraints. J Guidance Control Dyn 2008;31(3):543-53.

4. Jorris TR, Cobb RG. Three-dimensional trajectory optimization satisfying waypoint and no-fly zone constraints. J Guidance Control Dyn 2009;32(2):551-72.

5. Zhao J, Zhou R. Reentry trajectory optimization for hypersonic vehicle satisfying complex constraints. Chin J Aeronaut 2013;26(6): 1544-53.

6. Bollino KP. High-fidelity real-time trajectory optimizaiton for reusable launch vehicles [dissertation]. Monterey: Naval Postgraduate School; 2006.

7. Qi G, Kang W, Ross IM. A pseudospectral method for the optimal control of constrained feedback linearizable systems. IEEE Trans Autom Control 2006;51(7):1115-29.

8. Zhao J, Zhou R, Jin X. Progress in reentry trajectory planning for hypersonic vehicle. J Syst Eng Elect 2014;25(4):627-39.

9. D'Souza SN, Sarigul-Klijn N. Survey of planetary entry guidance algorithms. Prog Aerosp Sci 2014;68:64-74.

10. Xie Y, Liu LH, Tang GJ, Zheng W. Highly constrained entry trajectory generation. Acta Astronaut 2013;88(7-8):44-60.

11. Wingrove RC. Survey of atmosphere re-entry guidance and control methods. AIAA J 1963;1(9):2019-29.

12. Zimmerman C, Dukeman G, Hanson J. Automated method to compute orbital reentry trajectories with heating constraints. J Guidance Control Dyn 2003;26(4):523-9.

13. Xue S, Lu P. Constrained predictor-corrector entry guidance. J Guidance Control Dyn 2010;33(4):1273-81.

14. Joshi A, Sivan K, Amma SS. Predictor-corrector reentry guidance algorithm with path constraints for atmospheric entry vehicles. J Guidance Control Dyn 2007;30(5):1307-18.

15. Xu ML, Chen KJ, Liu LH, Tang GJ. Quasi-equilibrium glide adaptive guidance for hypersonic vehicles. Sci China Technol Sci 2012;55(3):856-66.

16. Lu P, Forbes S, Baldwin M. Gliding guidance of high L/D hypersonic vehicles. Reston: AIAA; 2013, Report No.: AIAA-2013-4648.

17. Harpold JC, Graves Jr CA. Shuttle entry guidance. J Astronaut Sci 1979;27(3):239-68.

18. Harpold JC, Gavert DE. Space shuttle entry guidance performance results. J Guidance Control Dyn 1983;6(6):442-7.

19. Mease KD, Chen DT, Teufel P, Schonenberger H. Reduced-order entry trajectory planning for acceleration guidance. J Guidance Control Dyn 2002;25(2):257-66.

20. Saraf A, Leavitt JA, Chen DT, Mease KD. Design and evaluation of an acceleration guidance algorithm for entry. J Spacecraft Rockets 2004;41(6):986-96.

21. Leavitt JA, Mease KD. Feasible trajectory generation for atmospheric entry guidance. J Guidance Control Dyn 2007;30(2):473-81.

22. Shen Z, Lu P. Onboard generation of three-dimensional constrained entry trajectories. J Guidance Control Dyn 2003;26(1): 111-21.

23. Shen Z, Lu P. Dynamic lateral entry guidance logic. J Guidance Control Dyn 2004;27(6):949-59.

24. Li HF, Zhang R, Li ZY, Zhang R. New method to enforce inequality constraints of entry trajectory. J Guidance Control Dyn 2012;35(5):1662-7.

25. Mease KD, Kremer JP. Shuttle entry guidance revisited using nonlinear geometric methods. J Guidance Control Dyn 1994;17(6): 1350-6.

26. Lu P. Entry guidance and trajectory control for reusable launch vehicle. J Guidance Control Dyn 1997;20(1):143-9.

27. Morio V, Cazaurang F, Vernis P. Flatness-based hypersonic reentry guidance of a lifting-body vehicle. Control Eng Pract 2009;17(5):588-96.

28. Dukeman G. Profile-following entry guidance using linear quadratic regulator theory. Reston: AIAA; 2002, Report No.: AIAA-2002-4457.

29. Lu P. Regulation about time-varying trajectories: precision entry guidance illustrated. J Guidance Control Dyn 1999;22(6):784-90.

30. Fahroo F, Ross I. Trajectory optimization by indirect spectral collocation methods. Reston: AIAA; 2000, Report No.: AIAA-2000-4028.

31. Tian B, Zong Q. Optimal guidance for reentry vehicles based on indirect Legendre pseudospectral method. Acta Astronaut 2011;68 (7-8):1176-84.

32. Park C, Guibout V, Scheeres DJ. Solving optimal continuous thrust rendezvous problems with generating functions. J Guidance Control Dyn 2006;29(2):321-31.

33. Peng HJ, Gao Q, Wu ZG, Zhong WX. Optimal guidance based on receding horizon control for low-thrust transfer to libration point orbits. Adv Space Res 2013;51(11):2093-111.

34. Han J. From PID to active disturbance rejection control. IEEE Trans Ind Elect 2009;56(3):900-6.

35. Vincent J, Morris D, Usher N, Gao ZQ, Zhao S, Nicoletti A, et al. On active disturbance rejection based control design for superconducting RF cavities. Nucl Instrum Methods Phys Res Sect A 2011;643(1):11-6.

36. Phillips T. A common aero vehicle (CAV) model, description and employment guide. Arlington: Schafer Corporation for AFRL and AFSPC; 2003.

Guo Jie received his Ph.D. degree from Beijing Institute of Technology

in 2010. His current research area is flight mechanics and control.