Mars atmospheric entry guidance for optimal terminal altitude

Accepted Manuscript Mars atmospheric entry guidance for optimal terminal altitude Jiateng Long, Ai Gao, Pingyuan Cui, Yang Liu PII:

S0094-5765(17)30905-0

DOI:

https://doi.org/10.1016/j.actaastro.2018.12.001

Reference:

AA 7223

To appear in:

Acta Astronautica

Received Date: 3 July 2017 Revised Date:

23 October 2018

Accepted Date: 1 December 2018

Please cite this article as: J. Long, A. Gao, P. Cui, Y. Liu, Mars atmospheric entry guidance for optimal terminal altitude, Acta Astronautica (2019), doi: https://doi.org/10.1016/j.actaastro.2018.12.001. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Mars Atmospheric Entry Guidance for Optimal Terminal Altitude Jiateng Longa, b, c, 1, Ai Gao a, b, c, 2, Pingyuan Cui a, b, c, 3* and Yang Liu a, b, c, 4 School of Aerospace Engineering, Beijing Institute of Technology, 100081 Beijing, China b Key Laboratory of Autonomous Navigation and Control for Deep Space Exploration, Ministry of Industry and Information Technology, 100081 Beijing, China c Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, 100081 Beijing, China 1 Ph.D. Candidate; [email protected]. 2 Associate Professor; [email protected]. 3 Corresponding author, Professor; [email protected]. 4 Post Doctor; [email protected].

RI PT

a

Abstract

SC

Maximizing the terminal altitude for Mars atmospheric entry has long been investigated on trajectory design to allow a sufficient timeline margin for subsequent operations and the scientific requirements of exploring Mars ancient highland. The purpose of this paper is to design an onboard Mars atmospheric entry guidance algorithm, which can achieve the optimal

M AN U

terminal altitude at the predetermined terminal flight range. Two important characters that the optimal bank angle profile are revealed in this paper, offering the gateway to the application of the optimal guidance by onboard parameter searching, which will be accomplished by the numerical predictor-corrector strategy. Moreover, suboptimal situations are also investigated considering the performance restriction of reactive control system (RCS). Effectiveness of the proposed guidance algorithm is demonstrated using scenarios of the Mars Science Laboratory (MSL) mission.

1. Introduction

TE D

Keywords: Mars atmospheric entry; terminal altitude; optimal guidance; numerical predictor-corrector strategy

Mars atmospheric entry guidance is to modulate a single variable, i.e., the bank angle, through which the vehicle can be decelerated appropriately for parachute deployment and delivered to an accurate location for the following descent and landing sequence with the assistance of aerodynamic force. Due to the highly thin atmosphere of the Mars (about only 1% in

EP

density to that of the Earth), compared with Earth reentry, in order to sufficiently decelerate and maneuver the vehicle, a relatively lower altitude (denser layer of atmosphere) of the vehicle needs to be achieved to generate enough aerodynamic forces. As a result, this usually lead to the difficulty in delivering the vehicle to a high elevation on Mars surface. Therefore,

AC C

one of the most challenging objectives for future Mars landing exploration mission is to deliver the vehicle to a higher terminal altitude at the predetermined landing site to meet the requirements of scientific interests, sufficient timeline margin and safety for the subsequent operations [1]. Because of such technological interest on achieving higher terminal altitude in Mars atmospheric entry, trajectory optimization is a beneficial and highly focused way for the offline atmospheric entry trajectory design. Researches in [2-8] studied the trajectory optimization in Mars atmospheric entry with respect to the optimal terminal altitude. In these literature, the characteristics of the optimal entry trajectory illustrated on the mission design parameters[2, 3], reachable footprint[5], controllable set [6], etc. Among the trajectory optimization methods, evolutionary heuristic algorithms like genetic algorithm (GA), particle swarm optimization (PSO) [2-4] are widely applied because of the ability of global optimal solution searching. However, the drawbacks on the huge computational load make it impossible to be adopted onboard.

ACCEPTED MANUSCRIPT On the other hand, in 2011, Mars Science Laboratory (MSL) mission, the first successful Mars mission to perform a guided entry [9], indicates a tendency of applying online closed-loop guidance in the future Mars atmospheric entry. In the design of MSL entry guidance, the lift of the vehicle is manipulated to balance the requirement of minimizing range error and ensuring a safe terminal altitude [9]. In fact, the Earth reentry of the Apollo capsules[10] and the space shuttle [11] have already provided beneficial knowledge for Mars atmospheric entry guidance algorithm design. From then on, tremendous

RI PT

efforts have been made to improve both categories of the Earth reentry guidance: the reference trajectory tracking method [12] (the two branches of which, i.e. trajectory planning [13] and trajectory tracking [14], were also researched in detail) and the predictive guidance method [15, 16]. The most recent Orion mission motivated the research of reentry guidance one step further and verified the method of the Apollo reentry guidance once again [16]. All of these works are highly beneficial for the guidance algorithm design of planetary entry including that of Mars atmospheric entry.

SC

With the experiences of Earth reentry, Carman et al. [17] applied the Apollo guidance to Mars precise landing, which is the baseline scheme of the MSL entry guidance. The two categories of Earth reentry guidance turn out to be still effective. Relevant contributions [17-19] also have been made on both of the reference trajectory tracking and predictive guidance

M AN U

algorithm for Mars atmospheric entry. In [20], Kluever investigated the performance of the these two categories of Mars entry guidance with respect to the landing accuracy. In addition, various studies have focused on Mars precision landing concerning on the specifics of Mars atmosphere. The algorithm in [21] dealt with high model uncertainties of the Mars predictive entry guidance. Most these efforts have been on the landing accuracy, which leave a potential to enhance the Mars atmospheric entry guidance algorithm by introducing onboard optimal guidance with respect to the terminal altitude. In fact, onboard optimal guidance have been introduced to different mission backgrounds [22-24]. Generally, the theoretical foundation of the onboard optimal guidance lies in the boundary value problem (BVP), which is derived from the

TE D

calculus of variations and the Pontryagin’s minimum principle (PMP). It is indicated by PMP that the fuel-optimal and timeoptimal control are of bang-bang structure. If similar conclusion can also be drawn in the Mars atmospheric entry with respect to the optimal terminal altitude as well, it will be of great benefit for both onboard guidance and the offline entry trajectory optimization to realize the maximal terminal altitude by reducing the computational load significantly. In this paper, the idea of optimal guidance is introduced to onboard guidance algorithm design for Mars atmospheric

EP

entry to explore the behavior of the optimal bank angle in view of the maximal terminal altitude. To this end, the form of the optimal entry guidance is firstly investigated by PMP with the constraints of terminal flight range, transforming the

AC C

optimality problem into an onboard parameter-searching problem. Suboptimal situations are also investigated considering the constraints of bank angle rate and acceleration due to the performance restriction of reactive control system (RCS). The effectiveness of the purposed algorithm will be illustrated with the background of MSL atmospheric entry. The contribution of this paper are three-fold: 1) with the help of the onboard predictor-corrector strategy, the proposed optimal guidance is capable of achieving the optimal terminal altitude in Mars atmospheric entry and the terminal position accuracy as well. 2) theoretically derived from the PMP, two important characters of the optimal guidance command profile are revealed. These characters guarantee the optimality of the guidance command, offering the gateway to the onboard parameter searching of the proposed optimal guidance. Moreover, the offline trajectory optimization will also benefit from these theoretical results by transforming the optimization to the switching instant searching, given the difficulties of the initial guess of the costate variables with indirect optimization method and the computational load with the direct method. 3)

ACCEPTED MANUSCRIPT with the optimality analysis, the proposed optimal guidance algorithm does not rely on any offline trajectory optimization, which saves the onboard storage and also shows the flexibility under the situation of in-flight uncertain environment.

2. Problems in Mars atmospheric entry In Mars atmospheric entry phase, the motion of the entry vehicle is governed by the entry dynamics with respect to a stationary atmosphere of a rotating planet. Entry dynamics is described by the following dimensionless equations [25]:

θɺ =

v cos γ sin ψ r cos φ

φɺ =

v cos γ cosψ r

RI PT

rɺ = v sin γ

SC

vɺ = -D − g sin γ + Ω2 r cos φ ( sin γ cos φ − cos γ sin φ cosψ )

 Ω 2 r cos φ  v2  − g  cos γ  + ( cos γ cos φ − sin γ sin φ cosψ ) v  r  

1 v

M AN U

γɺ =  L cos σ + 

 Ω2r 1  L sin σ v 2 + cos γ tan φ sin ψ  + sin φ cos φ sin ψ − 2Ω ( tan γ cos φ cosψ − sin φ ) v  cos γ r  v cos γ

ψɺ = 

(1) (2) (3) (4) (5)

(6)

where r is the Mars’ geocentric distance, normalized by the radius of the planet R0 , θ and φ are the longitude and latitude respectively, ψ is the heading angle measured from north clockwise, v is the planet-relative velocity normalized by

vc = g0 R0 with g0 = 3.71 m s2 being the gravitational acceleration at Martian surface, γ is the flight-path angle. σ is

TE D

2 the bank angle. Ω is the rotating angular velocity of Mars. g = µ r is the local the gravitational acceleration normalized

by g0 with µ being the gravitational constant of the Mars. The drag and lift acceleration D and L are D=

q 1 2 S ρv CD = d , β 2 m

L=

1 2 S ρv CL = D ⋅ ( L D ) 2 m

(7)

EP

Both of the drag and lift accelerations are normalized by g0 , where CD and CL are drag and lift coefficient, respectively, S 2 is the reference area of the vehicle, m is the vehicle mass, ρ is the Martian atmospheric density, qd = ρ v / 2 is the

AC C

dynamic pressure, β = m / SCD is the ballistic coefficient, and L/D is the lift-to-drag ratio of vehicle. The differential equations Eqs. (1)-(6) are with respect to the dimensionless time τ = t

R0 g0 .

For the Mars atmospheric density, a conventional exponential model is used

ρ = ρ0 e

h0 − h hs

(8)

−4 3 where ρ0 = 2 ×10 kg m is the reference density, h = r − R0 is the altitude of the entry vehicle, h0 = 40000 m is the

reference height, hs = 7500 m refers to the atmospheric scale height. The entry motion can be decoupled as longitudinal and lateral motion by ignoring the planet rotation ( Ω=0 ). Let s be the the range-to-go from the vehicle position to the predetermined parachute deployment position along the great arc of planet, and its derivative with respect to dimensionless time τ is

ACCEPTED MANUSCRIPT sɺ = − v cos γ

(9)

Then, by defining the state variable x = [ r , v,γ , s ] and combining Eqs. (1), (4), (5) and (9), the longitudinal model is T

v sin γ     -D − g sin γ    xɺ = f ( x, σ ) =  L cos σ  v g  +  −  cos γ   r v  v    −v cos γ

RI PT

shown in Eq. (10) ignoring the planet rotation

(10)

The terminal constraint of Mars atmospheric entry is that the flight range should reach the predetermined terminal range

s v(t ) =v − s f = 0 f f

SC

when the parachute is deployed, which can be written as

(11)

where t f is the free terminal time when parachute deployment condition is satisfied, s f is the predetermined terminal range

M AN U

(in general is zero), v f is the predetermined parachute deployment velocity.

3. Optimality of the entry guidance with respect to the terminal altitude With the longitudinal model Eq. (10) and terminal constraints Eq. (11) of the Mars atmospheric entry, the optimal terminal altitude problem is described as follows. Problem 1 (optimal terminal altitude problem):

According to PMP, the performance index is with the following form

TE D

f J = Φ  x ( t f ) , t f  + ∫ V  x ( t ) , u ( t ) , t  dt t0

t

(12)

In order to maximize the terminal altitude, the Mayer type of the performance index is introduced as J = −r ( t f

)

(13)

EP

where, Φ  x ( t f ) , t f  = − r ( t f ) , V  x ( t ) , u ( t ) , t  = 0 . With the longitudinal model Eq. (10) and the performance index Eq.(13), the Hamiltonian can be obtained as

AC C

L  v g H = λ T xɺ = λr v sin γ + λv (-D − g sin γ ) + λγ  cos σ +  −  cos γ  + λs ( −v cos γ ) r v  v 

(14)

T where λ =  λr , λv , λγ , λs  is the vector of costate variables, which are governed by Euler-Lagrange condition λɺ = − ∂H ∂x ,

i.e.,

 D 2 g sin γ + r  hs

  L cos σ  v 2 g  cos γ  + −  + λγ    r v  r    vhs 2λ D  L cos σ  1 g   λɺv = −λr sin γ + v − λγ  +  + 2  cos γ  + λs cos γ 2 v v r v    

λɺr = −λv 

λɺγ λɺs

v g = −λr v cos γ + λv g cos γ + λγ  −  sin γ − λs v sin γ r v  =0

(15)

ACCEPTED MANUSCRIPT The boundary conditions are shown in Eq. (16), and transversality conditions in Eqs. (17)-(18), respectively.

v0

γ0

v ( t f ) − v f  T =0 s0 ] , G  x ( t f ) , t f  =  s (t f ) − s f    0  −1  0  υ1  υ1  = 0  υ2   0     1 υ2 

 −1 0  0  1 ∂Φ ∂G + λ (t f ) = υ=  +  0  0 ∂x ( t f ) ∂x ( t f )     0  0 T

 ∂Φ ∂G T H (t f ) = −  +  ∂t ∂t f  f

(16)

(17)

RI PT

x ( t0 ) = x0 = [ r0

 ν =0   t =t f

(18)

structure

( cos σ )

 ( cos σ )min , if λγ > 0 = ( cos σ ) max , if λγ < 0

or

∗

 σ

σ =

 σ

, if λγ > 0 , if λγ < 0 min

max

M AN U

∗

SC

In order to minimize the Hamiltonian Eq. (14), the bank angle σ , which is the guidance command, should be the bang-bang

(19)

where the λγ is the switching function. Noted that cos σ will be undetermined when the switching function λγ ≡ 0 , which will lead to a singular optimal control problem. For this problem, detailed proof in the Appendix A shows that the singular optimal control is not possible for the bang-bang structure in Eq. (19). Based on such facts, the bang-bang structure of Eq. (19) is shown in Fig. 1, where the switching vectors t seven = ts(1)

ts( 2 ) ⋯ ts( 2 n )  and t sodd = ts(1) T

ts( 2) ⋯ ts( 2 n −1)  are for T

vector will be denoted as t s = ts(1)

TE D

the situations of even and odd number of switches, respectively. In the following discussion, for generality, the switching

ts( 2) ⋯ ts( N −1)

ts( N )  , where N is the number of switches. T

It should be noted that Mars atmospheric entry is terminated when the parachute deployment condition (i.e., the predetermined deployment velocity v f ) is reached. Therefore, the terminal time t f is free and the searching interval

EP

t0 , t f  of the switching vector t s is unfixed. Considering onboard parameter searching, a parameter with fixed searching

interval will be beneficial. To this end, relative velocity v can be selected as the searching parameter because of the fixed

AC C

interval i.e.,  v f , v0  , when the boundary conditions are given. Therefore, the switching velocity v s is chosen as the switching vector for the following discussion.

Fig. 1 Illustration of Eq. (19) with the conditions of even and odd number of switches

In the bang-bang structure Eq. (19), the switching function λγ , which is a costate variable, is critical for determining of the optimal bank angle profile. However, it is well known that even for offline trajectory optimization, solving the costate variables of BVP has long been a problem due to the difficulty in initial value guess. This makes the calculation of the costate practically impossible for onboard purpose. In fact, according to Eq. (19), the optimal bank angle σ

∗

can be fully determined by the following two characters of λγ :

ACCEPTED MANUSCRIPT 1) The sign of λγ . The sign of λγ determines which bound the optimal bank angle σ

∗

should take, i.e., the problem of

switching sequence. This problem will be analyzed in Section 3.1. 2) The roots of λγ . The roots of λγ , which is denoted as the switching vector, determine when optimal bank angle σ

∗

should switch, i.e. the problem of switching instant. This problem will be discussed in Section 3.2.

RI PT

Based on these two characters of λγ , when the switching sequence is determined, the problem of onboard optimal guidance can be converted into the problem of onboard searching of the switching vector v s = vs(1)

vs( 2) ⋯ vs( n )  . In the T

following discussion, the trajectory with respect to optimal terminal altitude can be denoted as

(

xɺ * = f x* ; σ *

)

is the optimal bank angle profile.

(20)

SC

where x* is the optimal entry trajectory, σ

*

M AN U

3.1 Determination of switching sequence

To fully determine the switching sequence, the number of switches for the optimal terminal altitude problem is essential. To this end, the longitudinal model Eq. (10) is equivalently transformed into the linear form zɺ = Az + Bu [26], where 0 0 A= 0  0

z2

z3

z4 ] = [ r T

rɺ s sɺ ] = [ r v sin γ

r   sin γ  u   ɺɺ u =  1 =   =  ɺɺ u  2   s   − cos γ

0 0  0  1

(21)

s −v cos γ ] . The control vector is

T

TE D

The state vector is z = [ z1

1 0 0 0 1 0 0 0  ,B =  0 0 0 1   0 0 0 0

v cos γ   vɺ   sin γ = v sin γ  γɺ   − cos γ

T

fv ( z )  v cos γ     v sin γ   f1 ( z ) u0 + f 2 ( z ) 

(22)

EP

 f ( z ) sin γ + f 2 ( z ) v cos γ   cos γ  = v +  f1 ( z ) vu0  − f v ( z ) cos γ + f 2 ( z ) v sin γ   sin γ 

where f v ( z ) = vɺ = − D − g sin γ , f1 ( z ) = L v , f 2 ( z ) = ( v r − g v ) cos γ , and u0 = cos σ . It should be noted that f1 ( z ) > 0

AC C

along the trajectory. Then, the velocity v and flight-path angle γ can be obtained as

v = rɺ 2 + sɺ2 = z 22 + z42 tan γ = −

(23)

z rɺ =− 2 sɺ z4

Accordingly, Problem 1 can be described as Problem 1′ .

Problem 1′ (optimal terminal altitude problem): The performance index is

J = − r ( t f ) = − z1 ( t f The Hamiltonian is

)

(24)

ACCEPTED MANUSCRIPT H z = λˆ T zɺ = λ1 z2 + λ2 u1 + λ3 z4 + λ4 u2

( )

( )

(25)

= H 0 λˆ, z + H1 λˆ, z u0 where

( ) H ( λˆ, z ) = ( λ cos γ + λ sin γ ) f ( z ) v = ( λ

H 0 λˆ, z = λ1 z2 + λ3 z4 + ( λ2 sin γ − λ4 cos γ ) f v ( z ) + ( λ2 cos γ + λ4 sin γ ) f 2 ( z ) v 2

4

1

The boundary conditions are

z ( t0 ) = z0 = [ z10

z20

z30

z40 ] = [ r0 T

v0 sin γ 0

s0

and transversality conditions are

T

 =0  

0 −1    z2 ( t f )υ z1   0  2 2  υ z1   z2 ( t f ) + z4 ( t f =   υz2 1  υ z 2     z4 ( t f )υ z1  0   z2 t + z2 t   2 ( f ) 4 ( f

M AN U

0   z2 ( t f )   −1   0   z22 ( t f ) + z 42 ( t f ∂Φ z ∂GzT λˆ ( t f ) = + υz =   +  0  0 ∂z ( t f ) ∂z ( t f )    z4 ( t f ) 0   z2 t + z2 t  2 ( f ) 4 ( f

−v0 cos γ 0 ]

SC

v ( t f ) − v f   z 22 ( t f ) + z42 ( t f ) − v f =   Gz  z ( t f ) , t f  =   s (t f ) − s f   z3 ( t f ) − s f   

RI PT

1

(26)

2 + λ4 tan γ ) f1 ( z ) v cos γ

) )

 ∂Φ z ∂GzT  H z (t f ) = −  + νz  =0  ∂t  ∂t f  f  t =t f

(28)

(29)

TE D

where λˆ = [ λ1

    )     ) 

(27)

ɺ T λ2 λ3 λ4 ] is the costate vector, which is governed by Euler-Lagrange conditions λˆ = − ∂H z ∂ z shown

in Eq. (30)

EP

λɺ1 = 0, λɺ2 = −λ1 , λɺ3 = 0, λɺ4 = −λ3

(30)

where [ λ10

AC C

Therefore, the solutions of Eq. (30) are

λ1 ( t ) = λ10 , λ2 ( t ) = λ20 − λ10 t

(31)

λ3 ( t ) = λ30 , λ4 ( t ) = λ40 − λ30 t

T λ20 λ30 λ40 ] is the initial value of λˆ .

According to the Eqs. (14) , (17) , (25) and (28), it can be concluded that

λr = λ1 , λv = λ2 sin γ − λ4 cos γ

(32)

λs = λ3 , λγ = ( λ2 cos γ + λ4 sin γ ) v

In Eq. (17) and Eq. (28),

λv ( t f ) = υ1 , λs ( t f ) = υ 2

λ2 ( t f ) =

z 2 ( t f ) υ z1

z 22 ( t f ) + z 42 ( t f

)

, λ3 ( t f ) = υ z 2 , λ4 ( t f ) =

z 4 ( t f )υ z 1

z 22 ( t f ) + z 42 ( t f

(33)

)

ACCEPTED MANUSCRIPT With Eq. (32) and Eq. (33), it can be derived that υ1 = υ z1 ,υ2 = υ z 2 . According to Eq. (30) and Eq. (32) that λɺr = λɺ1 = 0   D 2 g sin γ    D 2 g sin γ  along the trajectory. In Eq. (15), λɺr ( t f ) =  −υ1  + = 0 . Since  +   is not a constant along the h r r  s   t = t f  hs   trajectory, it can be concluded that υ1 = υ z1 = 0 . As a result, λ2 ( t f ) = λ4 ( t f ) = 0 .  ( cos σ )min , if H1 ( λˆ, z ) > 0 ∗ u0* = ( cos σ ) =  ˆ ( cos σ )max , if H1 ( λ, z ) < 0

or

∗

 σ

σ =

 σ

( )

RI PT

To minimize the Hamiltonian, the optimal bank angle should be of the bang-bang structure max min

, if H1 ( λˆ, z ) > 0

, if H1 ( λˆ, z ) < 0

(34)

Therefore, the sign of the switching function H 1 λˆ, z determines the optimal bank angle profile. In Eq. (26), denote

( )

( )

( )

SC

Hˆ 1 λˆ, z = λ2 + λ4 tan γ . Then, We obtain that H 1 λˆ, z = Hˆ 1 λˆ, z f1 ( z ) v cos γ . Since f1 ( z ) v cos γ = L cos γ > 0 along the

( )

trajectory regardless of the bank angle profile, the number of roots of Hˆ 1 λˆ, z also determines the number of switches

( )

M AN U

during entry. Moreover, according to Eq. (26) and Eq. (32), we obtain that H1 λˆ, z = λγ f1 ( z ) , indicating that the roots of

( )

Hˆ 1 λˆ, z is the roots of λγ , which also shows the equivalence of Problem 1 and 1′ . By substituting Eq. (28) into Eq. (31),

λ10 = −1 and λ30 = υz 2 , i.e. λ1 = −1 , λ2 = λ20 + t , λ4 = λ40 − υ z 2 t . Since λ2 ( t f ) = 0 , it can be derived that λ20 = −t f . Thus, the switching function

Hˆ 1 = λ2 + λ4 tan γ = ( λ20 + t ) + ( λ40 − υz 2t ) tan γ

TE D

(35)

Fig. 2 Illustration of Eq. (34) with the conditions of even and odd number of switches Proposition 1: For the bang-bang structure in Problem 1 and Problem 1′ , the minimum number of switches N = 2 .

EP

Proof: According to Eq. (34), the switch occurs when the switching function Hˆ 1 = 0 . Denote the switching instants as

t = ts(i ) ( i = 1, 2,⋯, N ) and ts(i ) < ts( j ) when i < j . Correspondingly, the tangent value of flight-path angle at each of the () witching instant are tan γ s , ( i = 1, 2,⋯, N ) . Then, the following relationship can be established at each switching instant

AC C

i

according to the switching function Eq. (35)

(

) (

)

− λ20 + t s(i ) = λ40 − υ z 2 t s(i ) tan γ s(i ) , ( i = 1, 2,⋯ , N )

(

( )

)

(36)

(

The right hand side of Eq. (36) indicates that each pair of ts( ) , λ4 ts( ) tan γ s( ) , i = 1, 2,⋯ , N along with t f , λ4 ( t f ) tan γ f i

i

i

)

locates on the line of −λ2 = − ( λ20 + t ) , whose slope kλ2 = −1 , i.e.,

( )

λ4 ( t f ) tan γ f − λ4 ts(i ) tan γ s(i ) (i )

t f − ts Since λ4 ( t f ) = 0 , Eq. (37) can be further written as

= −1, ( i = 1, 2,⋯ , N )

(37)

ACCEPTED MANUSCRIPT

(

)

− λ40 − υ z 2 t s( ) tan γ s( ) i

i

t f − t s( ) i

= −1, ( i = 1, 2,⋯ , N )

(38)

The relationship between optimal bank angle and λ2 , λ4 tan γ is illustrated in Fig. 2. Therefore, the solution of

λ  W  40  = b υ z 2  where

 ts(1) − t f   (2)   ts − t f ts( 2) tan γ s( 2 )    ⋮ ⋮   i () i i t s( ) tan γ s( )  , b =  ts − t f   ⋮ ⋮     N 1 − ( ) N −1 N −1 −tf ts( ) tan γ s( )  t s  (N ) (N ) (N )  ts tan γ s  N ×2  ts − t f

            N ×1

SC

ts(1) tan γ s(1)

M AN U

 − tan γ s(1)  ( 2)  − tan γ s  ⋮   W = − tan γ s(i )  ⋮   ( N −1)  − tan γ s  (N)  − tan γ s

RI PT

λ4 ( t ) = λ40 − υ z 2t is determined by the following linear equations (39)

(40)

( )

The condition which ensures the existence of the solution to Eq. (39) is that rank (W ) = rank Wɶ . For ease of discussion, denote Wɶ = W b  and rank ( • ) is the rank of matrix. The situations of the solution of Eq. (39) is analyzed as follows.

( )

1) When rank (W ) = rank Wɶ = 1 , Eq. (39) has infinitely many solutions. Possibilities can be divided into N = 1 and

TE D

N ≥ 2.

If N = 1 , there is only one switch during entry. Since Eq. (39) has infinitely many solutions, when the line

λ4 ( t ) = λ40 − υ z 2t changes, according to the switching function Eq. (35), among the terminal time t f , the profile of flight

EP

path angle γ and the switching instant t s (i.e., the roots of Hˆ 1 = 0 ), at least one of these three parameters will necessarily change. This indicates that there should be infinitely many optimal trajectories. Note that these optimal trajectories satisfy the terminal constraint Eq. (11). If the terminal altitude is not taken into account, it is proved in the Appendix A that the

AC C

trajectories satisfying the terminal flight range constraint Eq. (11) are finite. Therefore, when the terminal altitude is considered, the optimal trajectories that satisfy the terminal constraint Eq. (11) are also finite, which is contradictory to that there are infinitely many optimal trajectories. If N ≥ 2 , W and b can be written as  − tan γ 1 t1 tan γ 1  W  W = =  1 ,b =  W   N ×2 W  N ×2

t1 − t f   b  =   N ×1

b1  b    N ×1

(41)

1×2 N −1 × 2 N −1 ×1 where W1 ∈ ℝ , W ∈ ℝ ( ) , and b ∈ ℝ ( ) . Denoting that Wɶ 1 = W1 b1  , Wɶ = W b  , where the ith

( i = 1, 2,⋯, N −1)

(

ɶ row of the matrix W

)

i i i Wɶ ( ) = k ( )Wɶ 1 , k ( ) ∈ ℝ , i.e.,

( )

ɶ i is W ( ) . Since rank (W ) = rank Wɶ = 1 , it can be concluded that

ACCEPTED MANUSCRIPT tan γ 1

(i)

tan γ s As a result,

t1

t s( ) i

=

t1 tan γ 1

(i )

(i )

ts tan γ s

=

t1 − t f ts( i ) − t f

= k ( i ) , ( i = 1, 2,⋯ , N − 1)

(42)

= 1 , indicating that there is only one switch, which is contradictory to that N ≥ 2 .

( )

( )

situation that rank (W ) = rank Wɶ = 2 .

( )

RI PT

To summarize, both of the situations for rank (W ) = rank Wɶ = 1 are impossible, and the possibility only lies in the

2) When rank (W ) = rank Wɶ = 2 , Eq. (39) has a unique solution, and N ≥ 2 . W and b can be written as   ≜   N ×1

 b1  b   2  b  N ×1

(43)

SC

 − tan γ 1 t1 tan γ 1  W1   t1 − t f     W =  − tan γ 2 t2 tan γ 2  ≜ W 2  , b = t 2 − t f   N × 2  W  N × 2  b W

1×2 N − 2 ×2 N − 2 ×1 where W1 ,W2 ∈ℝ , W ∈ ℝ ( ) , b ∈ ℝ ( ) , W1 and W2 are linearly independent, b1 = t1 − t f , b2 = t2 − t f . Denoting

M AN U

ɶ ɶ i Wɶ1 = W1 b1  , Wɶ 2 = W2 b2  , Wɶ = W b  , where the ith ( i = 1, 2,⋯, N − 2 ) row of the matrix W is W ( ) . If N = 2 , i.e. W = 0 and b = 0 . The unique solution of Eq. (39) is

λ40 =

t2 b1 tan γ 2 − t1b2 tan γ 1 ( t1 − t2 ) tan γ 1 tan γ 2

υz 2 =

(44)

b1 tan γ 2 − b2 tan γ 1 ( t1 − t2 ) tan γ 1 tan γ 2

TE D

It can be verified that with the solution Eq. (44), the switching condition Hˆ 1 ( t1 ) = Hˆ 1 ( t2 ) = 0 is satisfied. Therefore, both of the line λ4 ( t ) = λ40 − υ z 2 t and optimal trajectory has a unique solution.

( )

If N > 2 , Eq. (39) still has the unique solution Eq. (44). Since rank (W ) = rank Wɶ = 2 , it can be concluded that

(

)

AC C

EP

i i i i i Wɶ ( ) = k1( )Wɶ1 + k 2( )Wɶ 2 , k1( ) , k 2( ) ∈ ℝ , i.e.,

 tan γ s(i ) = k1(i ) tan γ 1 + k2(i ) tan γ 2  ( i ) (i) (i) (i ) ts tan γ s = k1 t1 tan γ 1 + k 2 t2 tan γ 2  (i) (i) (i) ts − t f = k1 b1 + k2 b2

( )

( )

(45)

( )

i i i i In order to satisfy the switching condition Hˆ 1 ts( ) = λ2 ts( ) + λ4 ts( ) tan γ s( ) = 0 , the following condition should be

satisfied

bs(i )  ( b1 − b2 ) tan γ 1 tan γ 2 + ( b2 tan γ 1 − b1 tan γ 2 ) tan γ s(i )  + ( tan γ 2 − tan γ 1 ) b1b2 tan γ s(i ) = 0

(46)

() () i i where bs = ts − t f . For any specified N > 2 , according to Eq. (46), there are infinitely many combinations of k1( ) , k2( ) . i

i

Different combinations of k1( ) , k2( ) will result in different ts( ) (according to Eq. (45)) and different trajectories. Thus, for i

i

i

N > 2 , there is a unique solution of λ4 ( t ) = λ40 − υ z 2 t and there are infinitely many optimal trajectories. Therefore, it can be concluded that the minimum number of switches N = 2 .□

ACCEPTED MANUSCRIPT Remark 1: Although a lot of literatures on offline trajectory optimization draw the conclusion that the bank angle profile with respect to optimal terminal altitude is of the bang-bang structure, to the best knowledge of authors, none of them shows the details of such a bang-bang structure. In this Proposition, the number of switches of the optimal bank angle profile is demonstrated theoretically. In this sense, the conclusion of this Proposition will be beneficial to the offline trajectory optimization as well. Moreover, such a conclusion is also the gateway to the onboard application of this optimal guidance

RI PT

because the number of switches significantly influences the optimality of the bank angle profile. Surely, the guidance algorithm with one switch is the easiest to be conducted by the onboard predictor-corrector algorithm (which will be presented in detail in the next subsection). But as a cost, in this situation, the performance of terminal altitude will be severely undermined.

Remark 2: For Mars atmospheric entry, it should be noticed that frequent switching will prematurely deplete the propellant of the reaction control system. From this point of view, for onboard guidance algorithm, the number of switching N should

SC

take the minimum value which is necessary, i.e., N = 2 .

Remark 3: The conclusion that the minimum number of switches N = 2 is based on the constraints of both terminal

M AN U

velocity and terminal flight range. This conclusion may vary when the boundary conditions or performance index are different. For example, when only the terminal velocity is constrained and the terminal flight range is free, it can be verified that the number of switches N = 1 based on a similar way as illustrated above (where in this case, υ z 2 in Eq. (28) is zero). For this situation, same conclusion can also be referred in [3].

∗ Proposition 2: In the last switching interval t ∈ t s( 2 ) , t f  , the optimal bank angle σ should satisfy the following condition

σ

∗

2 t ∈t s( ) , t f   

=σ

(47)

min

TE D

Proof: In this proof, the longitudinal model in Problem 1 will be used. When t → t −f , since the λɺ ∈ C14×1 , with the terminal condition of the costate variable λ ( t f ) in Eq. (17), the limit of λɺ in Eq. (15)

AC C

EP

  D 2 g sin γ   L cos σ  v 2 g  cos γ   + −  −λv  +  + λγ   D 2 g sin γ       r  r v  r    −υ1  +  hs   vhs    r  hs        −λ sin γ + 2λv D − λγ  L cos σ +  1 + g  cos γ  + λs cos γ   2υ D 2 lim− λɺ = lim−  r  =  sin γ + 1 + υ2 cos γ  v v2 r v    t →t f t →t f  v    v g    υ v + g cos ( 1 ) γ − υ2 v sin γ  −λr v cos γ + λv g cos γ + λγ  −  cos γ − λs v sin γ     r v  0  t = t f    0  

(48)

Since υ1 = 0 , according to Eqs. (14) and (17) H ( t f ) =  −v sin γ − υ1 ( D + g sin γ ) + υ 2 ( −v cos γ )  t =t = − rɺ ( t f ) + υ 2 sɺ ( t f f

)

(49)

With the knowledge in Eq. (18) that H ( t f ) = 0 and the fact that rɺ ( t f ) < 0 and sɺ ( t f ) < 0 (for low lift entry vehicle,

γ ( t f ) < 0 ), the conclusion can be drawn that υ 2 > 0 along the trajectory. Then,

limλɺγ = ( v cos γ − υ2 v sin γ )t =t t →t −f

f

>0

(50)

ACCEPTED MANUSCRIPT Eq. (50) indicates that λɺγ > 0 is true in a certain left neighborhood of t = t f . Since λγ ( t f ) = 0 , we obtain that λγ < 0 in a certain left neighborhood of t = t f . □ With Proposition 2 and the nonexistence of singular optimal control, it can be concluded that for the optimal bank angle *

( t0 ) , the bound that it should take is related to the number of switches N and

σ

*

 σ

min

 σ

max

( t0 ) = 

, N is an even number , otherwise

Furthermore, with Proposition 1, the switching sequence with respect to time σ

sq

(51)

is minimum-maximum-minimum. The

SC

switching sequence with respect to time and relative velocity are shown in Fig. 3.

*

RI PT

in Fig. 1 and Fig. 2 at the initial time, i.e., σ

Fig. 3 Switching sequence with respect to time and relative velocity

M AN U

Remark 4: The offline trajectory optimizations in Ref. [2] and [3], in which the terminal range constraints are not considered, show that the optimal bank angle always tends to take the minimal value in the period before the parachute deployment in order to achieve the maximal terminal altitude. Such a tendency is because the entry vehicle can obtain the maximum component of the lift in the longitudinal plane to achieve the maximal terminal altitude, which is presented in Fig. 4. Eq. (47) shows that such tendency is still effective for optimal bank angle when considering the terminal condition requirement.

TE D

Fig. 4 Aerodynamic forces acted on the entry vehicle

3.2 Onboard switching velocity searching

With the problem of the switching sequence fully determined, the onboard guidance algorithm will focus on the problem of the switching instant, which means the switching velocity is to be searched onboard. The optimal bank angle profile σ

*

EP

with bang-bang structure Eq. (19) should satisfy the terminal constraints Eq. (11). For ease of the following discussion, the

( ) and v( ) = v ( t ( ) ) , which are to be searched by numerical predictor-

two switching velocities are denoted as vs(1) = v t s(1) corrector strategy [25], where ts( ) < ts(

2)

AC C

1

2

s

2

s

and vs( ) > vs( ) . 1

2

To ensure the onboard convergence robustness of the numerical predictor-corrector strategy, searching for one single guidance parameter is the common practice in state-of-the-art guidance algorithms [27]. It has been illustrated in Proposition 1 that there should be two switches to achieve the optimal terminal altitude, meaning that two variables are needed to be searched onboard. If treated as the multiple searching problem, it is likely to lead to a problem that only a local minimum is searched. To avoid such problem, it can be found that there are two terminal conditions needed to be satisfied, which are the terminal parachute deployment condition Eq. (11) and the maximal terminal altitude Eq. (13). In fact, vs( ) and vs( 1

separate for the problem in question: for a given vs( ) , the second switching velocity vs( 1

2)

2)

are not

is implicitly determined which can

be searched to satisfy the terminal parachute deployment condition Eq. (11). Therefore, the two variable searching problem

ACCEPTED MANUSCRIPT can be separated into a sequential searching problem with single parameter to be searched at one time, which is shown as follows: 1) The outer loop. In the outer loop, vs( ) is to be searched to satisfy the condition of optimal terminal altitude Eq. (13). 1

( )

( )

y1 vs(1) = −r vs(1) 1

2)

→ min

(52)

is searched to satisfy the condition of terminal parachute deployment

RI PT

2) The inner loop. For a given vs( ) , the vs(

v =v f

condition Eq. (11), i.e.

(

) ( )

y2 vs( 2) vs(1) = s vs( 2)

(

)

y2 vs( 2 ) vs(1) =

( ) ) = s ( v( ) )

where, y21 vs(

2

2

s

respect to vs( ) and vs( 1

2)

v =v f

v =v f

( )

(54)

respectively, and can be solved by a step-size-controlled Gauss–Newton method method illustrated

(1)  ( k )

=  vs 

( 2)  ( k +1)

vs  

( 2) ( k )

= vs 

(1)

−ζk

TE D

 vs  

2)

2

1 2 ( 2)  − s f  = y21 vs → min 2 

− s f . Therefore both of Eq. (52) and Eq. (54) can be treated as the univariate function with

(1)  ( k +1)

1

( )

1  (2) s vs 2 

in [27].

where, ζ k( ) and ζ k(

(53)

can also be searched by minimizing the following equation [25]

SC

2)

− sf = 0

M AN U

or alternatively, vs(

v =v f

( 2)

−ζ k

y1  vs(1) 

∂y1 vs(1) 

(k )

(k )

2 ∂y2 vs( ) 

(k )

(55)

∂vs(1)

∂vs(

2)

∂y v( 2)  ( k ) ∂v( 2)  s   21  s  

2

(56)

are the step-size parameters. Therefore, by onboard integrating Eq. (10) or Eq. (21) to predict the

closed-loop.

EP

terminal error of flight range, and solving Eq. (52) and Eq. (54) to correct this error, the optimal guidance in question is

AC C

4. Application of the optimal guidance to MSL atmospheric entry The proposed optimal guidance algorithm is further illustrate through numerical simulations with the background of entry phase of 2011 MSL mission. The boundary conditions are shown in Table 1 and Table 2 [6, 9, 28], respectively. In Table 1, the dispersion of all the initial conditions are normally distributed. The gravitational parameter and radius of Mars 4 3 2 are µ = 4.2828 ×10 km s and R0 = 3396 km , respectively. For the model of entry vehicle, the lift-to-drag ratio and the 2 ballistic coefficient are L D = 0.24 and β = 146kg m , respectively[28]. The uncertainties for both of them are set as

uniformly distributed with 10%.

Table 1 Nominal values and dispersions of initial conditions in simulations

ACCEPTED MANUSCRIPT

Entry altitude Longitude Latitude Velocity Flight path angle Heading angle

Value

Unit

3σ

125 126.73 -3.9186 6083.3 -15.4892 93.2065

km deg deg m/s deg deg

2.5 0.13 0.11 2.0 0.60 0.15

Dimensionless value 1.0371 -1.5721 -0.7662 1.6459 -0.2700 1.4421

RI PT

Parameter

Table 2 Terminal conditions in simulations

Longitude Latitude Velocity Flight range-to-go

Value

Unit

137.431 -4.552 403 0

deg deg m/s km

Dimensionless value 2.3981 -0.0794 0.1135 0

SC

Parameter

The nominal atmospheric density model in Eq. (8) and the dispersed atmospheric density model [29] are uitilized. Fig. 5

M AN U

demonstrates the 1000 runs of the dispersed atmospheric density model with respect to the altitude.

Fig. 5 1000 runs of the dispersed atmospheric density model with respect to the altitude

4.1 Nominal and Monte Carlo results

Nominal results of onboard optimal guidance are listed in Table 3 and Fig. 6 with different lower bound with the upper

TE D

bound being set as 70  . For comparison, the situation of single switch (i.e., the maximum-minimum switching sequence) is also listed as a baseline. vn is denoted as the velocity when the peak load factor is reached.

Table 3 Nominal results

σ

vs vn Peak load factor (m/s) (m/s) (deg) 10 (Baseline) 4225.9 16.551 4018.7 10 [1550.3, 3977.3] 13.703 4332.9 15 [1654.8, 3954.2] 13.776 4327.8 20 [1820.9, 3939.9] 13.878 4320.6 25 [2041.0, 3912.1] 14.010 4298.1 The situation of single switch shows the limited capability of achieving a higher terminal altitude. In contrast, for double

Terminal altitude (km) 11.478 12.202 12.187 12.139 12.027

AC C

EP

min

switches, the smaller the lower bound is selected, the higher terminal altitude can be reached, which can also be regard as an extension of the conclusion of Proposition 2.

Fig. 6 Nominal results of different bank angle profies

Fig. 7 Nominal results of states with respect to the relative velocity with different bank angle profiles

ACCEPTED MANUSCRIPT Fig. 7 shows the nominal results of the states with respect to the relative velocity with different bank angle profiles. It can be observed that with different bank angle profiles, the trajectory is capable of achieving the desired terminal position (see the figure of longitude, latitude and range-to-go for details). In order to verify the proposed algorithm under the situation of uncertain environments, the proposed optimal guidance algorithm is compared with the algorithm applied in MSL mission [17] and the recently proposed Gaussian Process Regression-based guidance [30] through Monte Carlo

RI PT

simulations. The initial states dispersion and the model uncertainties are considered. These two algorithms are chosen to make a comparison with ours is because the algorithm proposed in MSL mission is the one that is verified by the engineering application, while the Gaussian Process Regression-based guidance is a state-of-the-art near optimal guidance. In all the simulations, the guidance command is updated with frequency of 1Hz.

In the MC simulations, the boundary conditions in Table 1 and Table 2 are used for MC simulation. Fig. 8 illustrates the

SC

MC results of all of these three algorithms. It can be seen that similar terminal altitude (12.160 km in terms of the mean) are achieved by both of the algorithms proposed in this paper and that with the Gaussian Process Regression-based guidance. The altitude of the algorithm presented in [17], however, is about 1.2 km lower than that achieved by the previous two

M AN U

algorithms, demonstrating the effectiveness of the proposed guidance algorithm in achieving an optimal terminal altitude. As for the onboard computational load, all the simulations are carried out on the PC running a 64-bit Windows-7 operating system with a 3.80GHz processor and 4GB of RAM. The average prediction time of the kriging model in Gaussian Process Regression-based guidance is 0.763 ms, while for the proposed optimal guidance algorithm, such a number is 75.403 ms. It should be note that the Gaussian Process Regression-based guidance is based on the Kriging spatial statistical approach, which is essentially a guidance algorithm that interpolating a set of trajectories obtained by offline optimization according to the certain flight state. In order to cover the situations that vehicle may experience in the real flight situation

TE D

and show robustness, the database of these trajectories should be considerably large. This inherently requires a number of online memories that store these optimal trajectories and reduces the onboard computational burden. Instead of a set of optimal trajectories, the optimality of the proposed guidance algorithm are guaranteed by the conclusion derived from the PMP as Proposition 1 and 2. Therefore, there is no need for offline trajectory optimization in the proposed optimal guidance. As a result, the onboard computational load is relatively large but acceptable. Such a character is inherited from the state-of-

EP

the-art onboard numerical predictor-corrector strategy, which is also characterized with the ability of dealing with the

AC C

onboard off-nominal situations.

Fig. 8 Terminal altitude of the dispersion trajectories

4.2 Analysis of bang-bang structure by offline optimization In this subsection, the effect of different number of switches on terminal altitude are further investigated. To this end, single, triple and quintuple switches are conducted for the situation of odd number of switches, while double, quadruple and sextuple switches are for that of the even number of switches. For the offline optimization, the performance index Eq. (13), the boundary conditions Eq. (16) and parameters in Table 1 and Table 2 are utilized. The results of different number of switches are shown in Table 4. It can be observed that single switch is not optimal in terms of the terminal altitude maximization. As for the situation with double to sextuple switches, approximately the same terminal altitude can be achieved. Such a results confirms the conclusion of the Proposition 1 that the bang-bang structures

ACCEPTED MANUSCRIPT with the number of switches N ≥ 2 are capable of achieving the similar terminal altitude, and the bang-bang structure with double switches is the one with the minimum switches. The optimal bank angle profiles are shown in Fig. 9 in the category of odd and even number of switches.

Table 4 Result comparison under the bang-bang structure

5 6

Switching velocity (m/s) 4316.8 [1755.9, 3943.2] [2026.7, 4023.0, 6066.7] [2347.8, 3929.2, 4133.5, 4632.3] [2037.0, 4041.6, 5063.1, 5141.6, 6064.8] [1813.2, 2438.7, 2554.4, 2605.1, 2635.9, 4257.3]

Peak load factor 16.479 13.875 14.029 13.940 14.029 14.092

RI PT

2 3 4

Terminal altitude (km) 11.299 12.119 12.118 12.171 12.118 12.161

SC

Number of switches 1

Fig. 9 Bank angle of odd and even number of switches

M AN U

4.3 Suboptimal solution with RCS constraints

In state-of-the-art technology, the attitude (including the bank angle) of the vehicle is activated by the reactive control system (RCS). However, due to the RCS performance restriction, the ideal bang-bang structure is unrealizable in real entry flight. The rate and acceleration of the bank angle is constrained by

σɺ ≤ σɺ max , σɺɺ ≤ σɺɺ max

(57)

Therefore, an alternative suboptimal guidance is developed in this subsection for this problem by introducing a reference

TE D

bank angle profile

(

σ = f c, vs ;  σ

min

,σ

max



)

(58)

where vs is the switch velocity; c > 0 is the factor related to the constraints of σɺ and σɺɺ . For such a profile, the following

σ ( vs ) =

σ

dv

max

= σ =σ

+σ

min

2

dσ

AC C

EP

two conditions

max

(59)

dσ dv

=0 σ =σ

min

should be met to approximate the bang-bang profile. Then, the constraints in Eq. (57) will be satisfied by adjusting the parameter c onboard by the Gauss-Newton method demonstrated in Section V, in which σɺ and σɺɺ can be approximated with

σɺ =

σɺɺ =

d2 σ dt 2

=

dσ

d σ dv d σ = vɺ = σ ′ vɺ dv d t dv

(60)

d 2 σ  dv 2 d σ d 2 v = σ ′′ vɺ 2 + σ ′ vɺɺ   + dv 2  dt  dv dt 2

(61)

dt

=

ACCEPTED MANUSCRIPT where σ ′ and σ ′′ are the first and second order derivatives of σ with respect to relative velocity, which can be solved by Eq. (58), while vɺ and vɺɺ are evaluated with onboard inertial measurement unit (IMU) and finite difference method. In this paper, the following reference profile

max

min

− +

σ

max

−σ

min

1 + exp  c ( v − vs ) vm  σ max − σ min

1 + exp  c ( v − vs ) vm 

,σ ′ >0

(62) , σ′<0

or simply written as max

+σ

min

2

( )

σ − sign σ ′  

will be applied and further discussed, where vm = ( v0 + v f

)

min

−σ

max

2

+

σ

max

−σ

  vm  

min

1 + exp  c ( v − vs )

(63)

SC

σ

σ (v) =

RI PT

 σ  σ (v) =   σ 

2 . v0 and v f are the initial and terminal velocity, respectively.

The bank angle profile Eq. (63) will approach the bang-bang structure when c → +∞ , and dv

σ ′′ = where k1 = k 2 =

vm ( σ

c max

−σ

min

)

.

M AN U

dσ σ′= = − k1 ( σ − σ

d2 σ dv 2

min

)( σ

= −k2 ( 2 σ − σ

min

−σ

−σ

max

max

)

(64)

)

(65)

TE D

Then, the onboard bank angle command will be generated by Eq. (63) instead of Eq. (19) when the σɺ and σɺɺ constraints are imposed. By onboard searching the additional parameter c as well as the switching vector v s , the requirements of the rate and acceleration of bank angle can be satisfied. Fig. 10 illustrate the situations with different bank angle rate constraints (i.e., 15deg/s, 20deg/s, 25deg/s and 30deg/s,

angle is  20

70  .

EP

respectively) under the condition of double switch. The details of Fig. 10 are shown in Table 5, where the bounds of bank

AC C

Fig. 10 Suboptimal bank angle profile under the constraints of σɺ with double switches

Among the situations above with different levels of bank angle rate constraint, the bank angle suffers the toughest bank angle rate and acceleration at the first switch. The details in Table 5 also indicate that when bank angle rate and acceleration are constrained more strictly, the lower terminal altitude will be achieved.

Table 5 Results with different bank angle rate constraints under the condition of double switches

Maximum rate of bank angle (deg/s) 15 20 25

Switch velocity (m/s)

c

[1364.8, 3919.3] [1412.6, 3929.6] [1488.9, 3975.1]

[31.19, 22.96] [42.70, 27.90] [52.05, 28.08]

Maximum acceleration of bank angle (deg/s2) 7.03 12.23 18.65

Terminal altitude (km) 11.896 11.960 12.003

ACCEPTED MANUSCRIPT 30

[1972.7, 4388.0]

[61.36, 31.75]

27.62

12.100

5. Conclusion An onboard optimal guidance algorithm with respect to the maximum terminal altitude is proposed. The bang-bang structure is proven to be optimal for bank angle profile by PMP with terminal flight range constraint. Then, two propositions are put forward and proved, which are critical for the switching sequence determination. The proposed optimal guidance

RI PT

algorithm features the following aspects: Firstly, based on the predictor-corrector algorithm, with the proposed optimal guidance, optimal terminal altitude in Mars atmospheric entry will be achieved while satisfying the terminal position constraint. Secondly, two characters of the optimal guidance command profile derived from the PMP are essential to guarantee the optimality of the guidance command. These results facilitate the onboard application of the optimal guidance and benefit the the offline trajectory optimization as well. Thirdly, the proposed optimal guidance algorithm satisfies the

SC

optimality condition without the need of offline trajectory optimization, which saves the onboard storage and also shows the flexibility under the situation of in-flight uncertain environment. Simulation results show the effectiveness of the proposed optimal guidance method. Offline optimization is also conducted to analyze the effect of the number of switches, the

M AN U

conclusion of which confirms that of the Proposition 1. Furthermore, the suboptimal solution with control restrictions is proposed to avoid the bank angle rate and acceleration saturation.

Appendix A: Proof of the nonexistence of the singular optimal control

Assuming that singular optimal control exists in a finite interval [t1 , t2 ] ⊂ t0 , t f  , i.e. λγ ( t ) ≡ 0 , when t ∈ [t1 , t2 ] . Thus, in the time interval [t1 , t2 ] , λɺγ ( t ) ≡ 0 . According to Eqs. (30) and (32), λɺr = λɺ1 = 0 . By substituting λγ ( t ) ≡ 0 , λɺγ ( t ) ≡ 0

TE D

 D 2 g sin γ   D 2 g sin γ  and λɺr = 0 into the first of Eq. (15), −λv  +  ≡ 0 . Since the term  +  is not a constant during entry, r r  hs   hs 

λv should be zero and λɺv = 0 . Using these results, the second and third equation of Eq. (15) can be written as −λr sin γ + λs cos γ = 0 and λr v cos γ + λs v sin γ = 0 , respectively. Because cos γ ≠ 0 during entry, it can be concluded that

EP

λr (1 + tan 2 γ ) = 0 . Since (1 + tan 2 γ ) > 0 , λr = 0 in the time interval [t1 , t2 ] , which is contradict to the fact that λr = −1

AC C

during entry. Therefore, there is no existence of singular optimal control in Problem 1 and 1′ .□

Appendix B: Proof of the finiteness of the trajectories satisfying the terminal flight range constraint with single switch In this proof, the single switching velocity is denoted as vs . Therefore, according to the longitudinal model Eq. (10) v sɺ v cos γ dv = ∫ dv v 0 D + g sin γ vɺ v γ v ɺɺ γɺɺ dv γɺ ( v ) = ∫ dv = − ∫ v0 v v0 D + g sin γ ɺ

s (v) = ∫

v

v0

where

(66)

ACCEPTED MANUSCRIPT  1   r

γɺɺ =  L  −

 rɺ  g 1  1 g  L   g 2 − 2  cos σ + g  2 −  cos γ  +  2 +  cos γ + 2 cos σ  vɺ hs v  v r v v r v      

(67)

  sin γ D sin γ    2 g 2 sin γ + Dg D − g sin γ  = − L  + 2− +  cos σ +   cos γ  2 v hs  v r   r   

The derivatives of the variable upper limit definite integrals in Eq. (66) are

Γɺ d ( v ) =

ds ( v )

dv ɺ dγ ( v ) dv

=

d v v cos γ v cos γ sɺ = dv = ∫ v 0 dv D + g sin γ D + g sin γ vɺ

RI PT

Sd ( v ) =

d v γɺɺ γɺɺ γɺɺ =− ∫ dv = − = D + g sin γ vɺ dv v0 D + g sin γ

(68)

A smooth function of the bank angle profile with single switch in Eq. (62) is introduced as

max

derivative

ds ( v f dvs

) =  ds ( v )

dv dγɺ   dv dγɺ ( v ) dvs

 S (v) L = d sin σ  Γɺ d ( v ) v  Sd ( v f Γɺ ( v d

f

) (σ )

−σ

  S ( v ) L d cos σ = d   Γɺ ( v ) v dv s v = v f  d

(σ

−σ

min

  S (v) L dσ  sin σ = − d    Γɺ ( v ) v dvs v = v v = v f  d f

) exp c ( v − v ) v {1 + exp c ( v − v ) v }

max

−σ

s

min

m

2

s

f

min

f

s

(69)

as a function with respect to the switching velocity vs , the

m

) exp c ( v − v ) {1 + exp c ( v − v ) v }

max

max

1 + exp c ( v − vs ) vm 

TE D

=

)

σ

M AN U

Therefore, when treating the terminal flight range s ( v f

−

SC

σ (v) = σ

s

2

m

vm  cv s vm

 cvs   vm   v = v f

 L sin σ    v v = v f

(70)

min

Assuming that there are infinitely many trajectories that satisfy the terminal flight range constraint with single switch, i.e., there is an interval [ vs1 , vs 2 ] ⊂ v0 , v f  that when vs ∈ [ vs1 , vs 2 ] , s ( v f ) ≡ 0 . According to Eq. (70), this ds ( v f

) ≡0

EP

necessarily requires that

d vs

when vs ∈ [ vs1 , vs 2 ] , indicating that Γɺ d ( v f ) → ∞ , i.e. γɺɺ v = v → ∞ . However, f

AC C

according to Eq. (67), the condition that γɺɺ v = v → ∞ is impossible to be reached since both of the v f and r ( v f f

)

are

not zero. Therefore, there is no such an interval [ vs1 , vs 2 ] ⊂ v0 , v f  that when vs ∈ [ vs1 , vs 2 ] , s f ≡ 0 . Consequently, the trajectories that satisfy the terminal flight range constraint with single switch are not infinitely many.□

Fig. 11 Bank angle profiles with single switching

Acknowledgements This work was supported in part by the National Basic Research Program of China (973 Program) 2012CB720000, the National Natural Science Foundation of China 61374216, 61304226, 61304248, 61603039, the China Postdoctoral Science

ACCEPTED MANUSCRIPT Foundation 2016M591087, SAST Foundation SAST2016036, and the Science and Technology Innovation Team of Beijing Institute of Technology.

References

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

RI PT

SC

[5]

M AN U

[4]

TE D

[3]

EP

[2]

R. D. Braun and R. M. Manning, "Mars exploration entry, descent, and landing challenges," Journal of Spacecraft and Rockets, vol. 44, pp. 310-323, 2007. J. M. Lafleur and C. J. Cerimele, "Mars entry bank profile design for terminal state optimization," Journal of Spacecraft and Rockets, vol. 48, pp. 1012-1024, 2011. G. L. Jacob, G. Neeler, and R. Ramanan, "Mars Entry Mission Bank Profile Optimization," Journal of Guidance, Control, and Dynamics, pp. 1-12, 2014. M. J. Grant and G. F. Mendeck, "Mars science laboratory entry optimization using particle swarm methodology," in AIAA Atmospheric Flight Mechanics Conference and Exhibit, 2007. J. Benito and K. D. Mease, "Reachable and controllable sets for planetary entry and landing," Journal of guidance, control, and dynamics, vol. 33, pp. 641-654, 2010. J. Long, A. Gao, and P. Cui, "Controllable set analysis for planetary landing under model uncertainties," Advances in Space Research, vol. 56, pp. 281-292, 2015. P. Cui, Z. Zhao, Z. Yu, and J. Dai, "Terminal altitude maximization for Mars entry considering uncertainties," Acta Astronautica, vol. 145, pp. 446-455, 2018. S. Li and Y. Peng, "Mars entry trajectory optimization using DOC and DCNLP," Advances in Space Research, vol. 47, pp. 440-452, 2011. G. F. Mendeck and L. C. Mcgrew, "Entry Guidance Design and Postflight Performance for 2011 Mars Science Laboratory Mission," Journal of Spacecraft & Rockets, vol. 51, pp. 1094-1105, 2014. C. A. Graves and J. C. Harpold, "Apollo experience report: mission planning for Apollo entry," 1972. J. C. Harpold and C. A. Graves, "Shuttle entry guidance," in American Astronautical Society, Anniversary Conference, 25th, Houston, Tex., Oct. 30-Nov. 2, 1978, 35 p., 1978, pp. 239-268. H. Yan and Y. He, "Drag-tracking guidance for entry vehicles without drag rate measurement," Aerospace Science & Technology, vol. 43, pp. 372-380, 2014. K. Mease, D. Chen, P. Teufel, H. Sch-ograve, and nenberger, "Reduced-order entry trajectory planning for acceleration guidance," Journal of Guidance, Control, and Dynamics, vol. 25, pp. 257-266, 2002. S. Bharadwaj, A. V. Rao, and K. D. Mease, "Entry trajectory tracking law via feedback linearization," Journal of Guidance, Control, and Dynamics, vol. 21, pp. 726-732, 1998. A. I. Kozynchenko, "Predictive guidance algorithms for maximal downrange maneuvrability with application to low-lift re-entry," Acta Astronautica, vol. 64, pp. 770-777, 2009. S. Bairstow and G. Barton, "Orion Reentry Guidance with Extended Range Capability Using PredGuid," in AIAA Guidance, Navigation and Control Conference and Exhibit, 2007. G. L. Carman, Ives, D.G., and Geller, D.K., "Apollo-derived Mars precision lander guidance," in Atmospheric Flight Mechanics Conference, 1998. P. Richard W., "Numerical Roll Reversal Predictor-Corrector Aerocapture and Precision Landing Guidance Algorithms for the Mars Surveyor Program 2001 Missions," vol. 28, pp. 3107–3122, 1998. J. Benito, "Nonlinear Predictive Controller for Drag Tracking in Entry Guidance," in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2008. C. Kluever, "Entry guidance performance for Mars precision landing," Journal of guidance, control, and dynamics, vol. 31, pp. 1537-1544, 2008. A. I. Kozynchenko, "Analysis of predictive entry guidance for a Mars lander under high model uncertainties," Acta Astronautica, vol. 68, pp. 121-132, 2011. P. Lu, C. J. Cerimele, M. A. Tigges, and D. A. Matz, "Optimal Aerocapture Guidance," Journal of Guidance Control & Dynamics, vol. 38, pp. 553-565, 2015. A. I. Kozynchenko, "Development of optimal and robust predictive guidance technique for Mars aerocapture," Aerospace Science & Technology, vol. 30, pp. 150-162, 2013. H. Hu, S. Zhu, and P. Cui, "Desensitized optimal trajectory for landing on small bodies with reduced landing error," Aerospace Science and Technology, vol. 48, pp. 178-185, 2016. N. X. Vinh, A. Busemann, and R. D. Culp, Hypersonic and planetary entry flight mechanics: The University of Michigan Press, 1980.

AC C

[1]

ACCEPTED MANUSCRIPT

RI PT

SC

[30]

M AN U

[29]

TE D

[28]

EP

[27]

M. Zerar, F. Cazaurang, and A. Zolghadri, "Coupled linear parameter varying and flatness-based approach for space re-entry vehicles guidance," Iet Control Theory & Applications, vol. 3, pp. 1081-1092, 2009. P. Lu, "Entry guidance: A unified method," Journal of Guidance, Control, and Dynamics, vol. 37, pp. 713-728, 2014. S. Dutta and R. D. Braun, "Statistical Entry, Descent, and Landing Performance Reconstruction of the Mars Science Laboratory," Journal of Spacecraft & Rockets, vol. 51, p. S3, 2014. I. M. Meginnis, Z. R. Putnam, I. G. Clark, R. D. Braun, and G. H. Barton, "Guided Entry Performance of Low Ballistic Coefficient Vehicles at Mars," Journal of Spacecraft & Rockets, vol. 50, pp. 1047-1059, 2013. P. Ghosh and B. A. Conway, "Near-optimal feedback strategies synthesized using a spatial statistical approach," Journal of Guidance, Control, and Dynamics, vol. 36, pp. 905-919, 2013.

AC C

[26]

OJ

t0

t

V* V

V min

ts

2 n 1  2n s

t

tf

EP

2

AC C

ts

TE

D

max

t0 ts1

t0

t

V* V

Odd number of switches

M AN US C

OJ

RI PT

Even number of switches

t

max

V min

t0 ts1 ts 2 ts3

ts

2 n  2

ts

2 n 1

tf

t

RI PT

40

30 20 10

1000 |σ| ˙ |σ| ˙ |σ| ˙ |σ| ˙

2000

0

|σ| ˙ |σ| ˙ |σ| ˙ |σ| ˙

< 15deg/s < 20deg/s < 25deg/s < 30deg/s

2000

1000

< 15deg/s < 20deg/s < 25deg/s < 30deg/s

4000

5000

6000

7000

3000

4000

5000

6000

7000

3000

4000

5000

6000

7000

TE

60 40 20 0 −20 −40 0

1000

|σ| ˙ |σ| ˙ |σ| ˙ |σ| ˙

EP

−10 0

3000

< 15deg/s < 20deg/s < 25deg/s < 30deg/s

AC C

Bank angle rate, (deg/s)

20 0

Bank anlge accel. (deg/s2)

M AN US C

60

D

Bank angle, (deg)

80

2000

Relative velocity, m/s

RI PT M AN US C

V V max

AC C

EP

TE

vf

D

V min

vs

v0

v

RI PT

Even number of switches

Odd number of switches V*

V

V

M AN US C

V* max

V min

max

V min

t0

t

O4 tan J

t0

t

TE

D

O4 tan J

ts

2

ts

2 n 1  2n s

t

AC C

t0 ts1

EP

O2

tf

O2

t

t0 t 1 t  2 t 3 t  2 n  2 t  2 n 1 t f s s s s s

t

RI PT M AN US C

V

V

V max

V max

ts

2

tf

EP

1

t

TE

ts

AC C

t0

V min

D

V min

vf

vs

2

vs 1

v0

v

AC C EP TE

D

M AN US C

RI PT

RI PT M AN US C

120

100

TE

D

60

EP

40

20

0 −100

AC C

Altitude, km

80

−50

0

50

100

Deviation of Mars atmospheric density, %

150

200

RI PT

80

M AN US C

60

20

D

0

TE

−20

EP

−40 −60 −80 0

AC C

Bank angle σ, deg

40

1000

2000

3000 4000 Relative velocity, m/s

|σ|min=10deg |σ|min=15deg |σ|min=20deg |σ|

=25deg

min

5000

6000

7000

120

136

40

126 0

8000

−5 −10 −15

2000 4000 6000 Relative velocity v, m/s

8000

−4.3

−4.5 −4.6 0

2000 4000 6000 Relative velocity v, m/s

8000

|σ|min=10deg 600

100

95

90

85 0

−4.2

800

EP

Heading angle ψ, deg

0

−4.1

−4.4

8000

TE

105

2000 4000 6000 Relative velocity v, m/s

AC C

Flight path anlge γ, deg

130

D

2000 4000 6000 Relative velocity v, m/s

5

−20 0

132

128

20 0 0

134

Latitude φ, deg

60

−4

Range−to−go s, km

80

−3.9

M AN US C

100

−3.8

RI PT

138

Longitude θ, deg

Altitude h, km

140

|σ|min=15deg |σ|min=20deg

400

|σ|min=25deg

200 0

2000 4000 6000 Relative velocity v, m/s

8000

−200 0

2000 4000 6000 Relative velocity v, m/s

8000

600

600

RI PT

600 Min:10.580

500

Max:11.152

Mean:12.127

400

0

200

EP

200

100

300

TE

300

10.6 10.8 11 11.2 Terminal altitude, km (algorithm proposed in [17])

100

Max:12.381

Mean:12.160

400

Std:0.074

D

Number of cases

Std:0.079

AC C

Number of cases

400

500

Number of cases

Mean:10.852

Min:11.961

Max:12.343

M AN US C

500

Min:11.934

Std:0.076

300

200

100

0 11.8

12 12.2 12.4 Terminal altitude, km (algorithm proposed in [30])

0 11.8

12 12.2 12.4 Terminal altitude, km (algorithm proposed in this work)

Even number of switches

M AN US C

Odd number of switches 100

100

Single switch Triple switches Quintuple switches

60

TE

D

40

20

0 0

80

60

2000

EP

Bank angle, deg

80

Double switches Quadruple switches Sextuple switches

4000

Relative velocity, m/s

6000

40

20

0 0

2000

4000

Relative velocity, m/s

6000

ACCEPTED MANUSCRIPT Highlights

AC C

EP

TE D

M AN U

SC

RI PT

1) Optimal guidance achieves terminal position accuracy and optimal terminal altitude 2) Two important characters that the optimal bank angle profile are revealed 3) The optimal guidance algorithm does not rely on offline trajectory optimization