Trajectory optimization for asteroid landing with two-phase free final time

Journal Pre-proofs Trajectory optimization for asteroid landing with two-phase free final time Yingying Zhang, Jiangchuan Huang, Hutao Cui PII: DOI: Reference:

S0273-1177(19)30839-7 https://doi.org/10.1016/j.asr.2019.11.031 JASR 14555

To appear in:

Advances in Space Research

Received Date: Accepted Date:

3 July 2019 21 November 2019

Please cite this article as: Zhang, Y., Huang, J., Cui, H., Trajectory optimization for asteroid landing with twophase free final time, Advances in Space Research (2019), doi: https://doi.org/10.1016/j.asr.2019.11.031

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.

© 2019 Published by Elsevier Ltd on behalf of COSPAR.

Trajectory optimization for asteroid landing with two-phase free final time Yingying Zhang∗, Jiangchuan Huang, Hutao Cui School of Astronautics, Harbin Institute of Technology, Harbin, Heilongjiang150001, China

Abstract This work develops an autonomous trajectory planning algorithm for 6-DOF asteroid landing. The trajectory planning problem is formulated as a nonconvex time-optimal optimization problem with two-phase free final time, while the cost is regularized by augmenting a fuel consumption penalty. The nonconvex optimization problem is solved in successive solution method, so successive convexification is used to convert the original nonconvex problem into a sequence of convex subproblems, where each subproblem is obtained by linearizing the nonconvex dynamics and state constraints and using the velocity increment to give a convex expression of the fuel consumption penalty in cost function. Specifically, in the linearization, we divide the flight time interval into two parts and normalize each part using a time dilation coefficient to solve the problem that both the final times for the two flight phases are unknown, so that the original free final time problem turns to a fixed-time problem by minimizing the sum of the two time dilation coefficients and fuel consumption penalty. Besides, trust regions and virtual control are used to increase robustness of the algorithm. A convergence analysis is presented which indicates the successive solution will recover the local optimality of the original problem. Then the validity of the proposed algorithm and effects of different factors on flight time and fuel consumption are examined by simulations of landing on an irregular asteroid. Keywords: time-optimal, two-phase, free final time, convex optimization, successive solution, convergence analysis

∗ Corresponding

author Email addresses: [email protected] (Yingying Zhang), [email protected] (Jiangchuan Huang), [email protected] (Hutao Cui)

Preprint submitted to Journal of LATEX Templates

November 29, 2019

1. Introduction Asteroid exploration is an important way to learn about the origin and evolution of the solar system and understand asteroids that pose a potential threat to the Earth (Song and Gong, 2019; Ferrari and Lavagna, 2018). The current science-driven asteroid exploration mission focuses on landing and sampling, so that scientists can clarify the origin of life by analyzing samples acquired from a primordial celestial body. Hayabusa spacecraft is the first probe in the world to sample and return from an asteroid, it performed touching-downs on the asteroid Itokawa on November 2005 and returned less than a milligram of sample on June 2010 (Kawaguchi et al., 2008). Then Hayabusa 2 spacecraft successfully collected samples from the surface layer of the near-Earth asteroid 162173 Ryugu on February 2019 (Crane, 2019). OSIRIS-REx spacecraft will fly a Touch-And-Go trajectory to the asteroid 101955 Bennu to obtain a regolith sample from the surface and return in 2023 (Berry et al., 2013; Lauretta et al., 2017). Future asteroid exploration missions will need a high degree of landing accuracy, that is, more advanced and reliable trajectory design algorithms are needed to guide the precise and safe landing for the spacecraft. There are kinds of trajectory design algorithms for asteroid descent and landing, including closed-loop control based on current state (Guo et al., 2013; Furfaro et al., 2013; Lan et al., 2014a; Yang et al., 2017a), tracking a predesigned trajectory (Xiangyu et al., 2004; Zexu et al., 2012; Lan et al., 2014b; Li et al., 2006), optimal trajectory planning using the direct method (Lunghi et al., 2015; Hu et al., 2016) and the indirect method (Yang and Baoyin, 2015; Yang et al., 2017c). However, the closed-loop control significantly simplifies the problem compared to the optimal trajectory planning where the control constraint is considered, while the direct and indirect optimal trajectory planning methods are trapped in slow calculation speed and poor convergence respectively. In order to improve the efficiency of trajectory optimization, the optimal trajectory design problem is transformed into the convex optimization problem (Yang et al., 2017c). The convex optimization is first applied to the 3-DOF Mars landing problem (Acikmese and Ploen, 2007; Blackmore et al., 2010), and the relaxation technique of introducing slack variables is used to do lossless convexification for the nonconvex thrust magnitude and thrust pointing constraints (Acikmese and Ploen, 2007; Blackmore et al., 2010; Carson et al., 2011; Blackmore et al., 2012). The final optimal trajectory is generated by using the successive solution method, and the safe-guarding technique of trust regions and virtual control is applied in (Mao et al., 2016, 2017; Szmuk et al., 2016; Szmuk and Acikmese, 2018) to improve algorithmic robustness. Yang (Yang et al., 2017b) proposed a method for rapid 2

generation of time-optimal asteroid landing trajectories via convex optimization. Pinson (Pinson and Lu, 2018) proposed a propellant optimal trajectory design method for landing on irregularly shaped asteroids by convex optimization, and discussed the establishment of asteroid gravitational model and the lossless convexification for the thrust constraint where the thrust magnitude has lower and upper bounds. But both of the above methods are 3-DOF (degree-of-freedom) asteroid landing trajectory planning algorithms. The convex optimization problem is a class of optimization problems which can be solved efficiently, there are publicly and commercially available solvers with many published methods (Pinson and Lu, 2018), and has achieved great success on landing problems. This paper will use the convex optimization to do further research on the 6-DOF asteroid landing problem. The asteroid landing problem is formulated as a two-point boundary value problem, from a hovering point over the asteroid like 1 km (Kawaguchi et al., 2008) to the landing site. If the flight time is given, the problem can be solved as a propellant/fuel optimal problem or a minimum landing error problem. If the flight time is unknown , it can be done as a time-optimal problem with free final time. However, unless artificially set, the flight time of the optimal trajectory from the initial state to the terminal state with adjustable thrust is unknown. To determine the optimal flight time, in (Yang et al., 2017b) they proposed a minimum landing error problem to serve as a gateway for solving the time-optimal control problem and searched for the minimum flight time through a combination of extrapolation and bisection methods, in (Pinson and Lu, 2018) the convex optimization problem of minimizing the propellant was set as an inner loop and the optimal flight time was given by using a an outer optimization loop. But neither of the above two methods gives the optimal flight time directly, while generates the trajectory at the same time. In (Szmuk and Acikmese, 2018) the free final time problem was solved by normalizing the flight time into [0, 1] with a time dilation coefficient, then the flight time was given in a direct way by minimizing the time dilation coefficient. However, asteroids are highly irregular in shape and rotate quickly, the flight path of the vehicle should be constrained to avoid collisions with the asteroid surface. Based on (Liao-McPherson et al., 2016; Reynolds and Mesbahi, 2017), the landing process can be divided into two phases, each phase corresponds to a different state constraint, while the flight time in each phase is unknown. Thus, in this paper, we consider the translational motion and rotational motion of the vehicle to establish a 6-DOF dynamic model, and improve the time dilation coefficient method to apply it to a time-optimal trajectory planning problem with two-phase free final time. We will solve the prob-

3

lem under the convex optimization framework and further discuss the convergence of the successive solution method. The paper is organized as follows: Sec. II presents the time-optimal trajectory planning for asteroid landing with two-phase free final time, giving the dynamics and all constraints. Sec. III presents convex optimization framework for time-optimal asteroid landing trajectory planning, including the convexification and discretization for the original trajectory planning problem. Sec. IV presents the successive solution for generating the time-optimal landing trajectory, and gives the convergence analysis to the successive solution method, Sec. V presents some simulation results, Sec. VI concludes this paper.

2. Time-optimal trajectory planning for asteroid landing with two-phase free-finaltime The asteroid landing process can be considered in two phases which differ by whether the landing site is tracked, the descent phase and final landing phase (Liao-McPherson et al., 2016; Reynolds and Mesbahi, 2017). In the descent phase, the vehicle only need to move from the initial state to somewhere over the landing site, the flight time is [0, tc ]. In the final landing phase, the vehicle will land on the landing site from the state that the descent phase ends, while the landing site can be tracked from the beginning state, the flight time is (tc , tf ]. Where tc is the conversion time of the two phases. Here, both the transition time tc and the final time tf are unknown, tc can be regarded as the final time of the descent phase, so the time-optimal trajectory planning for asteroid landing we develop is a two-phase free final time optimization problem. Below we give the dynamics, all the state and control constraints, then give the formulation of this problem. 2.1. The dynamics First we define two coordinate frames (as shown in Fig. 1) related to the motion of the vehicle. FL is the asteroid-fixed coordinate frame centered at the asteroid’s center-of-mass, its Z-axis points along the spin axis of the asteroid, its X-axis points along the asteroid’s maximum (minimum) inertial principal axis, the Y-axis constitutes the right-handed system. FS is the vehicle-fixed coordinate frame centered at the vehicle’s center-of-mass, its z-axis points along the vehicle’s maximum inertial principal axis, its x-axis points along the vehicle’s minimum inertial principal axis, the y-axis constitutes the right-handed system. 4

Figure 1: Two coordinate frames diagram

The vehicle is considered to be a rigid body and equivalently actuated by three reaction wheels and six identical thrusters rigidly mounted to the vehicle axes. The fuel is considered to deplete at a rate proportional to the thrust magnitude, the mass depletion can be given by m(t) ˙ =−

|TS,x (t)| + |TS,y (t)| + |TS,z (t)| kTS (t)k1 =− Isp · g0 Isp · g0

(1)

where Isp is the vacuum specific impulse, and g0 is the Earth’s standard gravity constant. Then we express the translational motion in FL coordinate frame and the rotational motion in FS coordinate frame, and assume the position of the center-of-mass is constant. The translational and rotational motions are given by r˙ L (t) = vL (t)

(2)

v˙ L (t) = −2ω e × vL − ω e × (ω e × rL ) +

1 L C (σ LS )TS + ∇U (r) m S

1 R(σ SL )ω SSL 4   1 = R(σ SL ) ω − CSL (σ SL )ω e 4

σ˙ SL (t) =

˙ Jω(t) = MS (t) − ω × (Jω)

5

(3)

(4)

(5)



  R(σ SL ) =   CSL (σ SL )

1 + σ12 − σ22 − σ32

2(σ1 σ2 − σ3 )

2(σ1 σ3 + σ2 )

2(σ1 σ2 + σ3 )

1 − σ12 + σ22 − σ32

2(σ2 σ3 − σ1 )

2(σ1 σ3 − σ2 )

2(σ2 σ3 + σ1 )

1 − σ12 − σ22 + σ32

×



×

4 1 − σ TSL σ SL [σ SL ] − 8 [σ SL ] = E3 − !
2 1 + σ TSL σ SL !




2



0

  × [σ SL ] =  σ3  −σ2

     −σ3 0 σ1

σ2



  −σ1   0

In (2-5), the thrust TS (t) ∈ R3 , the torque MS (t) ∈ R3 , the vehicle’s position rL (t) ∈ R3 , the velocity vL (t) ∈ R3 , the gravity ∇U (r) ∈ R3 . The asteroid’s rotating angular velocity ω e ∈ R3 and the angular velocity ω ∈ R3 are expressed in FL coordinates and FS coordinates respectively. J ∈ S3++ is the the inertia tensor, σ SL = [σ1 , σ2 , σ3 ] is the MRPs (Modified Rodrigues Parameters) vector from FL to FS , CL S (σ LS ) is the transformation matrix between FS and FL coordinates, S L −1 T CSL (σ SL ) is the inverse transformation of CL = CL S (σ LS ), CL (σ SL ) = CS (σ LS ) S (σ LS ) .

2.2. The constraint Since the asteroid is irregularly shaped and rotating, a main constraint in the trajectory planning is to avoid collisions with the asteroid surface. In the descent phase, the landing site may not be visible by the vehicle’s landing camera, so that keeping the vehicle flying outside an ellipsoid enveloping the asteroid (as shown in Fig. 2) can meet the requirement of collision avoidance (LiaoMcPherson et al., 2016; Reynolds and Mesbahi, 2017). The three semi-major axes of the ellipsoid are given by a, b, c, the vehicle position rL (t) is simplified written as r(t), then the ellipsoid constraint can be given by 

1 a2

  1 ≤ r(t)T  0  0

0 1 b2

0

0



  T 0  r(t) = r(t) Re r(t), 

∀ t ∈ [0, tc ]

(6)

1 c2

In the final landing phase, to prevent a collision during landing, a glide-slope constraint introduced. The constraint is to restrict the vehicle to fly inside a cone centered at the landing site rtf with the conical angle θ measured from the normal vector external to the landing site n (Reynolds and Mesbahi, 2017) (as shown in Fig. 3). Here, both rtf and n are expressed in FL coordinate frame, thus the constraint can be given by

6

Figure 2: The ellipsoid constraint

T

cosθ ≤

[r(t) − rtf ] n kr(t) − rtf k

∀ t ∈ (tc , tf ]

(7)

In addition, in the final landing phase, landing site tracking is required to achieve pinpoint landing, so there is a constraint on the attitude of the vehicle to ensure that the landing site is in the FOV (field-of-view) of the landing camera (as shown in Fig. 4). Give the angle of FOV β, the installation position of vision sensor ρS , and the LOS (line-of-sight) direction dS , the attitude constraint which is also called the FOV constraint can be given by T  −CSL (σ SL )(r(t) − rtf ) − ρS dS cosβ ≤ k − CSL (σ SL )(r(t) − rtf ) − ρS kkdS k

∀ t ∈ (tc , tf ]

(8)

Another concern is the constraint on control, we constrain the thrust and the torque to lie in the magnitude interval [0, Tmax ] and [0, Mmax ], respectively. The thrust vector TS (t) and the torque MS (t)is simplified written to T(t) and M(t), then the control constraint can be given by kT(t)k∞ ≤ Tmax

kM(t)k∞ ≤ Mmax

∀ t ∈ [0, tf ]

(9)

Moreover, we restrict the mass of the vehicle to be greater than the dry mass mdry using the following the constraint mdry ≤ m(t),

∀ t ∈ [0, tf ]

7

(10)

Figure 3: The glide-slope constraint

2.3. Time-optimal trajectory planning problem for asteroid landing The time-optimal trajectory planning problem is based on solving a time-optimal control problem from the initial vehicle state at current time t0 to the terminal state at final time tf , that is min tf . In order to control the fuel consumption of the vehicle while obtaining tf , the cost function t

can be regularized by augmenting a fuel consumption penalty. The fuel consumption is the mass variation m(tf ) − m(0) which can be expressed by the thrust T(t) in Eq. (1), thus the cost function can be given by min tf +  t, T

where α =

1 Isp ·g0 ,

tf

Z

αkT(t)k1 dt

(11)

0

 is the weight of fuel consumption  ∈ R+ .

The summary of the time-optimal trajectory planning problem for asteroid powered landing with two-phase free final time is provide in Problem 1. Problem 1. Time-optimal trajectory planning problem with two-phase free final time (the original problem) min tf +  t, T

Z

tf

αkT(t)k1 dt

0

subject to m(t) ˙ = −αkT(t)k1

r˙ (t) = v(t)

¨r(t) = −2ω e × r˙ − ω e × (ω e × r) +

8

1 L C (σ)T + ∇U (r) m S

Figure 4: The field-of-view constraint

˙ σ(t) =

  1 R(σ) ω − CLS (σ)ω e 4

kT(t)k∞ ≤ Tmax

kM(t)k∞ ≤ Mmax

1 ≤ r(t)T Re r(t) kr − rtf k ≤ kCSL (σ)(r(t) − rtf ) + ρS k +

nT (r − rtf ) cosθ

˙ Jω(t) = M − ω × (Jω) mdry ≤ m(t)

∀ t ∈ [0, tf ]

∀ t ∈ [0, tc ] ∀ t ∈ [tc , tf ]

  S dTS CL (σ)(r(t) − rtf ) + ρS ≤ 0 ∀ t ∈ [tc , tf ] cosβkdS k

(12) (13)

r(0) = r0 , v(0) = v0 , σ(0) = σ 0 , ω(0) = ω 0 , m(0) = minit

(14)

r(tf ) = rtf , v(tf ) = vtf , σ(tf ) = σ tf , ω(tf ) = ω tf , mdry ≤ m(tf )

(15)

where Eq. (14) and Eq. (15) are the initial state and terminal state of the vehicle respectively. In Problem 1, all subscripts about FS and FL coordinate frames are removed for abbreviation.

3. Convex optimization framework for time-optimal asteroid landing trajectory planning The convex optimization problem is a class of optimization problems, it can be solved quickly and reliably in polynomial time with global optimality (Boyd and Vandenberghe, 2004), and there 9

are efficient algorithms such as interior-point method (IPM) (Dueri et al., 2016) to solve the convex problem. So convex optimization is considered to be a prime candidate for onboard automated guidance and real-time applications. Since the nonlinear dynamics in Eqs. (2-5) is not affine, the state constraints and cost funtion in Eqs. (6, 8, 11) are nonconvex, developing convexification for the Problem 1 to solve it under the convex optimization framework is necessary. The means of convexification is mainly to linearize the dynamics and nonconvex constraints and obtain their approximate convex expressions. Then we discretize the convex problem to cast the continuous-time optimizaation problem into a finite-dimensional parameter optimization problem. 3.1. Linearization Before linearization, we first define the state vector x(t) ∈ R12 and control vector u(t) ∈ R6 ,  T  T where x(t) = r(t)T , v(t)T , ω(t)T , σ(t)T , u(t) = T(t)T , M(t)T , so the dynamics in Eqs. (2-5) can be expressed as a vector-valued function

 T dx(t) ˙ T , ω(t) ˙ T , σ(t) ˙ T = X [x(t), u(t)] = r˙ (t)T , v(t) dt

(16)

Eq. (16) is a general expression of the dynamics, we give the following form to distinguish the dynamics in the descent phase from the dynamics in the final landing phase. dx(t) dx(t) dx(t) =ε + (1 − ε) dt dt dt

(17)

where ε = 1 when t ∈ [0, tc ], and ε = 0 when t ∈ (tc , tf ]. Eq. (17) means to artificially divide the dynamics into two parts, but the dynamics are essentially the same in the two flight phases, except that the initial state and the terminal state are different. The two flight phases are continuous, the terminal state of the descent phase is the initial state of the final landing phase. Since both tc and tf are unknown, we would like to express Eq. (17) in terms of normalized trajectory time, τ1 =

t tc

∈ [0, 1], τ2 =

t−tc tf −tc

∈ (0, 1], apply the chain rule to the left side of Eq. (17),

thus we obtain dx(t) dx(t) dτ1 dx(t) dτ2 =ε + (1 − ε) dt dτ1 dt dτ2 dt

(18)

Denote λ1 as the time dilation between τ1 and t, λ2 as the time dilation coefficient between τ2

10

and t, which are λ1 ,

dt dτ1

λ2 ,

dt dτ2

(19)

So according to the definition of τ1 and τ2 , we can get tc = λ1 and tf = λ1 + λ2 . Denote τ = τ1 or τ = τ2 , the dynamics in Eq. (16) can be rewritten in terms of τ by ˙ )= x(τ

dx(τ ) = [ελ1 + (1 − ε)λ2 ] X [x(τ ), u(τ )] dτ

(20)

Then we linearize the time-normalized dynamics, and approximate Eq. (20) by the first-order n o ˆ (τ ), u ˆ (τ ), m(τ Taylor expansion near the augmented reference trajectory x ˆ ), λˆ1 , λˆ2 . The lin-

earized dynamics can be given by

˙ ) = Gλ (τ ) [ελ1 + (1 − ε)λ2 ] + Gx (τ )x(τ ) + Gu (τ )u(τ ) + Gc (τ ) x(τ

(21)

ˆ (τ )] ∈ R12 Gλ (τ ) = X [ˆ x(τ ), u h i ˆ 1 + (1 − ε)λ ˆ 2 ∂X (x, u) xˆ (τ ),ˆu(τ ) ∈ R12×12 Gx (τ ) = ελ ∂x h i ∂X (x, u) ˆ 1 + (1 − ε)λ ˆ2 xˆ (τ ),ˆu(τ ) ∈ R12×6 Gu (τ ) = ελ ∂u Gc (τ ) = −Gx (τ )ˆ x(τ ) − Gu (τ )ˆ u(τ ) ∈ R12

In the state constraints, the ellipsoid constraint and the FOV constraint are nonlinear and nonconvex, they can also be turned to be convex by first-order Taylor expansion approximation. For the glide-slope constraint, it can be rewritten in the form of second-order cone which is convex. So the convex state constraints are by given by ˆ (τ1 ) − 2ˆ ˆ (τ1 )T Pe x x(τ1 )Pe x(τ ) ≤ 0 1+x kCe x(τ2 ) − Ce x(1)k ≤ Ce = [E3 , 03×9 ] ∈ R3×12

nT Ce [x(τ2 ) − x(1)] cosθ

(22) (23)

Pe = Ce T Re Ce ∈ R12×12

0 0 ˆ (τ2 ) ≤ 0 q [ˆ x(τ2 )] x(τ2 ) + q [ˆ x(τ2 )] − q [ˆ x(τ2 )] x

11

(24)

q (x) = kf (x)k +

dTS f (x) ∈ R cosβkdS k

T

0

q (x) =

f (x) dTS f (x)0 + f (x)0 ∈ R1×12 kf (x)k cosβkdS k

f (x) = CSL (σ)Ce x − CSL (σ)Ce xtf + ρS ∈ R3   ∂ CSL (σ)Ce x 0 S f (x) = CL (σ)Ce + Cσ ∈ R3×12 ∂σ Cσ = [03×9 , E3 ] ∈ R3×12 where x(1) is the time-normalized terminal state at tf , Pe and Ce are coefficient matrix. In fact, the glide-slope constraint can be approximated by a set of linear half-plane constraints (Richards and How, 2002), but considering that the problem is solved by second-order cone programming, so there is no linear approximation to Eq. (23). After time normalization, the mass depletion in Eq. (1) and cost function in Eq. (11) turn to m(τ ˙ ) = −α [ελ1 + (1 − ε)λ2 ] kT(τ )k1 min λ1 + λ2 + λ1  λ, T

Z

1

αkT(τ1 )k1 dτ1 + λ2 

0

Z

(25)

1

αkT(τ2 )k1 dτ2

(26)

0

3.2. Discretization To transform the original continuous-time optimal control problem into a finite-dimensional optimization problem, we need to discrete the normalized time domain and impose the constrains at the temporal notes. Discretize the normalized time domain [0, 1] into K1 equal intervals with the increment

1 K1

for the descent phase, and discretize the normalized time domain (0, 1] into

K2 equal intervals with the increment

1 K2

for the final landing phase. The temporal nodes are

denoted by τ1,k1 and τ2,k2 for the two phases respectively, τ1,k1 = k2 K2 ,

k1 K1 ,

k1 = 0, 1, 2, ..., K1 , τ1,k2 =

k2 = 1, 2, ..., K2 . In order to combine the two flight phases, we denote the temporal nodes in

the whole flight phase as τk , k = 0, 1, ..., K1 , K1 + 1, ..., K1 + K2 . Since τk must be incremented by k, we obtain that τk = τ1,k1 and ε = 1 when k1 = k = 0, 1, ..., K1 , τk = τ1,K1 + τ1,k2 and ε = 0 when k2 = k − K1 , k = K1 + 1, ..., K1 + K2 . Conduct the zero-order-hold control in each discrete

12

normalized time interval, thus we have  T x(k) = x(τk ) = r(τk )T , v(τk )T , ω(τk )T , σ(τk )T  T m(k) = m(τk ) u(k) = u(τk ) = T(τk )T , M(τk )T

(27)

Keep the thrust constant for τ ∈ [τk , τk+1 ), that is u(τ ) = u(k),

∀τ ∈ [τk , τk+1 ),

k = 0 , ..., (K1 + K2 − 1)

(28)

Then in the k th interval τ ∈ [τk , τk+1 ), the discrete expressions for mass depletion and dynamics can be given by   λ1 λ2 αkT(k)k1 m(k + 1) = m(k) − ε + (1 − ε) K1 K2

(29)

¯ x (k)x(k) + G ¯ u (k)u(k) + G ¯ λ (k) [ελ1 + (1 − ε)λ2 ] + G ¯ c (k) x(k + 1) = G

(30)

¯ k+1 − τk ) = e(τk+1 −τk )Gx (k) ¯ x (k) = φ(τ G Z τk+1 ¯ k+1 − ξ)Gu (k) dξ ¯ u (k) = φ(τ G ¯ λ (k) = G ¯ c (k) = G

τk τk+1

Z

τ Z kτk+1

¯ k+1 − ξ)Gλ (k) dξ φ(τ ¯ k+1 − ξ)Gc (k) dξ φ(τ

τk

Denote x(K) is the terminal state at tf , where K = K1 + K2 , so the discrete state constraints can be given by ˆ (k)T Pe x ˆ (k) − 2ˆ 1+x x(k)Pe x(k) ≤ 0 kCe x(k) − Ce x(K)k ≤

nT Ce [x(k) − x(K)] cosθ

0 0 ˆ (k) ≤ 0 q [ˆ x(k)] x(k) + q [ˆ x(k)] − q [ˆ x(k)] x

13

for k = 0, ... , K1

(31)

for k = K1 + 1, ... , (K1 + K2 )

(32)

for k = K1 + 1, ... , (K1 + K2 )

(33)

The continuous-time cost function in Eq. (26) can be discretized into

min

λ1 , λ2 , T

λ1 + λ 2 + 

K 1 −1 X k=0

λ1 αkT(k)k1 +  K1

K1 +K X2 −1 k=K1

λ2 αkT(k)k1 K2

(34)

The above-mentioned method of dealing with two flight phases solves the problem of unknown conversion time. But in Eq. (34), since dilation coefficients λ1 , λ2 and thrust T(k) are unknown, that

λ1 K1 αkT(k)k1

λ2 λ1 λ2 K2 αkT(k)k1 are nonconvex terms. Considering that K1 = ∆t1 and K2 = ∆t2 T(k)∆t2 1 so T(k)∆t m(k) and m(k) are velocity increments over a period of time ∆t, where

and

are time increments,

∆t represents ∆t1 or ∆t2 . Denote the velocity increment as ∆v(k) which can be given by ∆v(k) = v(k + 1) + ω e × r(k + 1) − v(k) + ω e × r(k) =

(35)

ˆ LS (k)] CL S [σ ˆ [r(k)] ∆tˆ T(k)∆t + ∇U m(k) ˆ

So the convex expression of the fuel consumption penalty can be deduce as

T(k)∆t =

(

CSL

"

ˆ ˆ ˆ [r(k)] ε λ1 + (1 − ε) λ2 ˆ SL (k)] ∆v(k) − ∇U [σ K1 K2

#)

m(k) ˆ

(36)

Then the cost function in Eq. (34) turns into

min

λ1 , λ2 , T∆t

λ 1 + λ2 + 

K1 +K X2 −1

αkT(k)∆tk1

(37)

k=0

Remark 1. The mass depletion in Eq. (29) is also nonconvex but not included in the state vector x. It is under consideration that the vehicle carries enough fuel and the constraint mdry ≤ m(k) can always be met. The planning for mass is neglected, the {m(k)} ˆ in the reference trajectory can o n ˆ1, λ ˆ2 . ˆ λ be computed by Eq. (29) with T(k),

After convexification and discretization, Problem 1 is put under the convex optimization fram-

work, and the convex optimization problem represented as Problem 2 is described as below. Problem 2. Discrete and convex time-optimal optimization problem

min

λ1 , λ2 , T∆t

λ 1 + λ2 + 

K1 +K X2 −1 k=0

14

αkT(k)∆tk1

subject to   λ2 λ1 αkT(k)k1 + (1 − ε) m(k + 1) = m(k) − ε K1 K2 ¯ x (k)x(k) + G ¯ u (k)u(k) + G ¯ λ (k) [ελ1 + (1 − ε)λ2 ] + G ¯ c (k) x(k + 1) = G ˆ (k) − 2ˆ ˆ (k)T Pe x x(k)Pe x(k) ≤ 0 1+x kCe x(k) − Ce x(K)k ≤

nT Ce [x(k) − x(K)] cosθ

0 0 ˆ (k) ≤ 0 q [ˆ x(k)] x(k) + q [ˆ x(k)] − q [ˆ x(k)] x

mdry ≤ m(k)

kT(k)k∞ ≤ Tmax

for k = 0, ... , K1 for k = K1 + 1, ... , (K1 + K2 ) for k = K1 + 1, ... , (K1 + K2 ) kM(k)k∞ ≤ Mmax

 T  T T m(0) = minit , x(0) = rT0 , v0T , ω T0 , σ T0 , x(K) = rTtf , vtf , ω Ttf , σ Ttf

(38) (39)

It must be noted that the problem 2 is a convex optimization problem, its solution has global optimality. But the Problem 2 is only the approximation of the Problem 1 which is a nonconvex optimization problem, the solution of Problem 2 is not in general the global optimal solution for Problem 1, it depends on the initial reference trajectory chosen. However, due to the deterministic behavior and numerical efficiency of convex optimization, Problem 2 can be treated as a subproblem of successive convexification for the onboard autonomous guidance algorithm.

4. Successive solution for generating the time-optimal landing trajectory In this section the successive solution method is adopted to solve the Problem 1 and generate the time-optimal trajectory. In order to increase robustness of the algorithm, two important techniques, virtual control and trust regions are introduced to eliminate the artificial infeasibility and capture the nonlinearity of dynamics and nonconvex state constraints in Problem 1. Then we will discuss the convergence property of the successive solution. 4.1. Successive solution For the successive solution method, the optimal solution is given by linearizing the dynamics and constraints successively and then solving a sequence of convex constrained subproblems, that is, Problem 2 is applied as a subproblem of successive convexification, each subproblem can be

15

solved as a second-order cone programming. Denote the ith successive solution by superscript i,  and the trajectory related parameters of the solution are expressed as xi (k), ui (k), mi (k), λi1 , λi2 .

Let i − 1th solution be the reference trajectory, that is

ˆ1 ˆ (k) ui−1 (k) , u ˆ (k) mi−1 (k) , m(k) xi−1 (k) , x ˆ λ1i−1 , λ

ˆ2 λ2i−1 , λ

i = 1, 2, 3, ... (40)

The dynamics and constraints used to find the ith solution are expressed by i − 1th solution can be given by

 i  λ λi−1 mi (k + 1) = mi (k) − ε 1 + (1 − ε) 2 αkTi−1 (k)k1 K1 K2   ¯ c (k) ¯ x (k)xi (k) + G ¯ u (k)ui (k) + G ¯ λ (k) ελi + (1 − ε)λi + G xi (k + 1) = G 1 2 1 + xi−1 (k)T Pe xi−1 (k) − 2xi−1 (k)Pe xi (k) ≤ 0

for k = 0, ... , K1

(41) (42) (43)

 nT Ce  i x (k) − xi (K) for k = K1 + 1, ... , (K1 + K2 ) (44) cosθ  0    0 q xi−1 (k) xi (k) + q xi−1 (k) − q xi−1 (k) xi−1 (k) ≤ 0 for k = K1 + 1, ... , (K1 + K2 ) (45) kCe xi (k) − Ce xi (K)k ≤

4.2. Trust regions and Virtual Control In order to solve Problem 1 under convex optimization framework, the nonlinear dynamics and nonconvex constraints are linearized to do convexification. This measure needs to be able to capture the nonlinearity of original problem, if not, it may make the problem unbounded. To eliminate the risk, we augment trust regions acting like penalty funtions into the cost function in Problem 2 to ensure the linearized trajectory does not deviate significantly from the previous iteration (Mao et al., 2016; Szmuk and Acikmese, 2018). The deviation between the ith and the i − 1th solution is δλi2 = λi2 − λi−1 2

δλi1 = λi1 − λi−1 1

δxik = xi (k) − xi−1 (k)

δuik = ui (k) − ui−1 (k) (46)

Then, define ∆iλ , ∆ix,k , ∆iu,k and impose the following constraint !

δλi1


2

!
2 + δλi2 ≤ ∆iλ

 i T i δxk δxk ≤ ∆ix,k

 i T i δuk δuk ≤ ∆iu,k

(47)

where ∆iλ , ∆ix,k , ∆iu,k are the square of the trust region radius respectively, the corresponding trust !
i−1 i−1 region centers are located at λi−1 , x (k) and ui−1 (k). 1 , λ2 16

T  Denote ∆i = ∆ix,0 , ... , ∆ix,K1 +K2 , ∆iu,0 , ... ∆iu,K1 +K2 −1 , ∆iλ , the cost function with trust

regions in the ith iteration can be given by

min ϕi [λ1 , λ2 , T∆t] := λi1 + λi2 + 

K1 +K X2 −1

i

αk (T(k)∆t) k1 + W∆ k∆i k1

(48)

k=0

Another concern is the linearization may not feasible, and the artificial infeasibility is often encountered during the first few the successive iterations (Szmuk and Acikmese, 2018) which can affect the iteration process and convergence (Mao et al., 2016). To eliminate this possibility, virtual controls Vxi (k) ∈ R12 , Vei (k) ∈ R+ and Vfi (k)+ ∈ R are added to the dynamics in Eq. (42) and constrains in Eqs. (43, 45) respectively, thus we can obtain   ¯ c (k) + Vi (k) ¯ x (k)xi (k) + G ¯ u (k)ui (k) + G ¯ λ (k) ελi + (1 − ε)λi + G xi (k + 1) = G x 1 2 Vei (k) ≥ 0

for k = 0, ... , K1

(50)

1 + xi−1 (k)T Pe xi−1 (k) − 2xi−1 (k)Pe xi (k) ≤ Vei (k) Vfi (k) ≥ 0

(49)

(51)

for k = K1 + 1, ... , (K1 + K2 )

(52)

 0    0 q xi−1 (k) xi (k) + q xi−1 (k) − q xi−1 (k) xi−1 (k) ≤ Vfi (k)

(53)

The penalty terms of virtual controls are added to the cost function, Eq. (48) turns to

min Li [λ1 , λ2 , T∆t] := λi1 + λi2 + 

K1 +K X2 −1

i

αk (T(k)∆t) k1

(54)

k=0

+ W∆ k∆i k1 + Wg

K1 +K X2 −1

kVxi (k)k1

k=0

+ Wh

"K 1 X

k=0

Vei (k) +

KX 1 +K2 K1 +1

#

Vfi (k)

where W∆ is the weight for the trust regions, Wg and Wh are penalty factors. The function of penalty terms is to assign a large objective function value to the point which is infeasible or attempts to cross the feasible region boundary, and force the minimal point of penalized problem to approach the feasible region infinitely or keep moving in the feasible region

17

until it converges to the minimal point of the original constrained optimization problem. It can be proved that as long as Wg ≥ 1 (Mao et al., 2016) and Wh ≥ 0 (Mao et al., 2018), the stationary point of problem described by Eqs. (49-54) is also the stationary point of Problem 2. In practical application, the values of Wg and Wh are very different, Wg should take a large value, while Wh takes a mall value. h iT i i Denote ∆mi = αk (T(0)∆t) k1 , ..., αk (T(K1 + K2 − 1)∆t) k1 as the fuel consumption      T T T matrix, λi = λi1 , λi2 , Vxi = Vxi (0)T , ... , Vxi (K1 + K2 − 1)T , Vei = Vei (0) , ... , Vei (K1 ) , h iT Vfi = Vfi (K1 + 1) , ... , Vfi (K1 + K2 ) , so the cost function in Eq. (54) can be rewritten as  min Li [λ, ∆m] := kλi k1 + k∆mi k1 + Wg kVxi k1 + Wh kVei k1 + kVfi k1 + W∆ k∆i k1

(55)

Finally, the convex discrete-time optimization subproblem for successive solution is provided in Problem 3. Problem 3. Convex finite-dimensional time-optimal optimization subproblem  min Li [λ, ∆m] := kλi k1 + k∆mi k1 + Wg kVxi k1 + Wh kVei k1 + kVfi k1 + W∆ k∆i k1 subject to  i−1  λi−1 λ αkTi−1 (k)k1 mi−1 (k + 1) = mi−1 (k) − ε 1 + (1 − ε) 2 K1 K2   ¯ x (k)xi (k) + G ¯ u (k)ui (k) + G ¯ λ (k) ελi + (1 − ε)λi + G ¯ c (k) + Vi (k) xi (k + 1) = G 1 2 x Vei (k) ≥ 0

for k = 0, ... , K1

1 + xi−1 (k)T Pe xi−1 (k) − 2xi−1 (k)Pe xi (k) ≤ Vei (k) kCe xi (k) − Ce xi (K)k ≤

 nT Ce  i x (k) − xi (K) cosθ

Vfi (k) ≥ 0

for k = K1 + 1, ... , (K1 + K2 )

for k = K1 + 1, ... , (K1 + K2 )

 0    0 q xi−1 (k) xi (k) + q xi−1 (k) − q xi−1 (k) xi−1 (k) ≤ Vfi (k)

δλi1 = λi1 − λi−1 1

δλi2 = λi2 − λi−1 2

δxik = xi (k) − xi−1 (k) δuik = ui (k) − ui−1 (k)

18

!

δλi1


2


2 ! + δλi2 ≤ ∆iλ

λ1,min ≤ λ1 ≤ λ1,max

mdry ≤ mi (k)

 i T i δuk δuk ≤ ∆iu,k

 i T i δxk δxk ≤ ∆ix,k

λ2,min ≤ λ2 ≤ λ2,max

kTi (k)k∞ ≤ Tmax

kMi (k)k∞ ≤ Mmax

 T  T T mi (0) = minit , xi (0) = rT0 , v0T , ω T0 , σ T0 , xi (K) = rTtf , vtf , ω Ttf , σ Ttf

(56) (57) (58)

where λ1,min and λ1,max , λ2,min and λa,max are constraints for time dilation coefficients λ1 and λ2 , which are used to limit the flight time in each flight phase for the vehicle. 4.3. Convergence Analysis   Define y = x(0)T , ..., x(K1 + K2 )T , u(0)T , ..., u(K1 + K2 − 1)T , λ1 , λ2 ∈ RN , where N =

12(K1 + K2 + 1) + 6(K1 + K2 ) + 2, and yi expresses the ith successive solution. Define hi (y) = T T   i i )T , where hi1,k [x(k)] = h1,0 , ..., hi1,K1 , hi2,K1 +1 , ..., hi2,K1 +K2 and gi (y) = (g0i )T , ..., (gK 1 +K2 −1  0     1+xi−1 (k)T Pe xi−1 (k)−2xi−1 (k)Pe xi (k), hi2,k [x(k)] = q xi−1 (k) xi (k) − xi−1 (k) +q xi−1 (k) ,   ¯ c (k). ¯ x (k)xi (k) − G ¯ u (k)ui (k) − G ¯ λ (k) ελi − (1 − ε)λi − G gki [x(k), u(k), λ1 , λ2 ] = xi (k + 1) − G 1 2

According to Eqs. (49-53), when hij (y) ≤ 0, j = 0, ..., K1 , K1 + 1, ..., K1 + K2 , it must have Vei (k) = 0 or Vfi (k) = 0, and when hij (y) > 0, it must have Vei (k) = hi1,k [x(k)] , k = 0, ..., K1 or

Vfi (k) = hi2,k [x(k)] , k = K1 + 1, ..., K1 + K2 . So the new expression for Eq. (55) is min Li (y) := kλi k1 + k∆mi k1 + Wg kgi (y)k1 + Wh

KX 1 +K2 j=0

 max 0, hij (y)

(59)

  Denote si (y) = si0 , ... , sij , ... , siK1 +K2 as other convex constraint vector, define a set F , where  F = y | sij (y) ≤ 0, j = 0, 1, ..., K1 + K2 .

According to Eqs. (6, 20, 24), discrete forms of the original nonconvex dynamics and state

constraints can be expressed as x(k+1) = X [x(k), u(k), λ1 , λ2 ], 1−x(k)Pe x(k) ≤ 0 and q [x(k)] ≤

19

0. If the dynamics and state constraints are not linearized, Eq. (59) should be the form as below.

min J i (y) := λi1 + λi2 + 

K1 +K X2 −1

i

αk (T(k)∆t) k1

(60)

k=0

+ Wg

K1 +K X2 −1 k=0

+ Wh

(K 1 X

k=0

i  

x (k + 1) − X xi (k), ui (k), λi1 , λi2 1

) KX 1 +K2    max 0, q xi (k) max 0, 1 − xi (k)Pe xi (k) + 

K1 +1

Define di = yi − yi−1 , so the trust region term W∆ k∆i k1 in Eq. (55) can be replaced by !
T W∆ di di . Since the first order Taylor expansion is applied to dynamics in Eq. (30), ellipsoid constraint in Eq. (31) and FOV constraint in Eq. (33), so that Eq. (60) can be expressed by di  T !
T Li (d) , Li−1 (y + d) := J i−1 (y) + ∇J i−1 (y) di + W∆ di di

(61)

Define the actual change of cost in Eq. (60) between the i − 1th iteration and ith iteration as ∆J i (y) = J i−1 (y) − J i (y) = J i−1 (y) − J i−1 (y + d)

(62)

Define the predicted change by linear approximation between the i − 1th iteration and ith iteration as  T !
T ∆Li (d) = J i−1 (y) − Li (d) = − ∇J i−1 (y) di − W∆ di di

(63)

Then the ratio of actual change and predicted change can be given by

 T !
T − ∇J i−1 (y) di − di ∇2 J i−1 (y + ξd)di ∆J i (y) γ = = T T ∆Li (d) − [∇J i−1 (y)] di − W∆ (di ) di i

(64)

Now we prove that an infinite successive solution sequence {yi } is generated by solving Problem ¯ which is also the stationary point of the Problem 1, it 3 successively, and {yi } has limit points y ¯ , and k∇J i (y)k → 0, γ i → 1. means when kdi k → 0 that yi → y ¯ . Use proof by contradiction, First, we prove that there is a subsequence yi that converges to y assume that when kdi k → 0, for all sufficiently large i, k∇J i (y)k ≥ ε, ε is a positive number, γ i does not converge to 1.

20

 0  T i−1 ¯ i = 0, we can get d ¯ i = − ∇J (y) , according to Eq. Let Li (d) = − ∇J i−1 (y) − 2W∆ d 2W∆   ¯ i is also the solution of ∆Li (d) 0 = 0, it means that the optimal solution of Eq. (61) (63), d

¯ i ∈ F , the maximum predicted change will make the maximum predicted change. When yi−1 + d is ∆Li (d) =

k∇J i−1 (y)k2 4W∆

dih ∈ F and ∆Li (d) ≥ 1 i−1 (y)kkdi k 2 k∇J

=

1 i−1 ¯ i k, (y)kkd 2 k∇J

1 i−1 (y)kkdih k. 2 k∇J

if not, there must a certain dih , so that yi−1 +

Based on the assumption, we can obtain ∆Li (d) ≥

≥ 12 εkdi k for all sufficiently large i, so for kγ i − 1k

W !di
T di − !di
T ∇2 J i−1 (y + ξd)di ∆

kγ i − 1k =

T i T i i−1 i

− [∇J (y)] d − W∆ (d ) d



!
T !
T

W∆ di di + di ∇2 J i−1 (y + ξd)di ≤ 1 i 2 εkd k



i 2 ! i
T 2 i−1

∇ J (y + ξd) di W∆ d + d ≤ 1 i 2 εkd k

o

n

2

!
T = W∆ di + 2 di ∇2 J i−1 (y + ξd) ε

(65)

Based on Eq. (65), when kdi k → 0, kγ i − 1k → 0, γ i → 1, that means the assumption is not ¯ and k∇J i (y)k → 0. valid, there must be a subsequence yi converges to y ¯ . Use proof by contradiction, Next we prove that there is no subsequence yi does not converge to y suppose that there are such a subsequence yi and a positive value ε, k∇J im (y)k ≥ ε for im . According to the previous proof, k∇J i (y)k has a subsequence that converges to 0. So there is an index set {in }, for each in there is k∇J in (y)k ≤ 13 ε. And for im ≤ i ≤ in , there are k∇J i (y)k ≥ 31 ε and ∆Li (d) ≥ 61 εkdi k. Combine Eqs. (62)-(64), we can obtain J i−1 (y) − J i (y) ≥

1 i i−1 εγ ky − yi k 6

(66)

So for kyim − yin k there is  1 i im 1 εγ ky − yin k ≤ εγ i kyim − yim+1 k + kyim+1 − yim+2 k + ... + kyin−1 − yin k (67) 6 6       ≤ J im (y) − J im+1 (y) + J im+1 (y) − J im+2 (y) + ... + J in−1 (y) − J in (y) ≤ J im (y) − J in (y)

21

When i → ∞, the left-hand of Eq. (67) tends to 0, so the right-hand also tends to 0. Therefore when kyim − yin k → 0, k∇J im (y) − ∇J in (y)k → 0. Take a large enough index i to ensure k∇J im (y) − ∇J in (y)k ≤ 13 ε, and we can deduce ε ≤ k∇J im (y)k = k∇J im (y) − ∇J in (y) + ∇J in (y)k

(68)

≤ k∇J im (y) − ∇J in (y)k + k∇J in (y)k ≤

1 1 2 ε+ ε= ε 3 3 3

The derivation is contradictory, the assumption is not valid. So the infinite sequence {yi } ¯ of the original generated by the successive solution method will converge to the stationary point y ¯ is not in general a global optimum, problem. It is important to point out that the stationary point y the convergence analysis can only indicate that the successive convexification method will recover local optimality for the original problem. Solving a convex subproblem in the successive solution process has good reliability and ensures the successive iteration can proceed steadily after the initial iteration starting, which is very important for the feasibility of onboard real-time applications. 5. Simulation and Results In this section, the simulation results are presented to demonstrate feasibility and validity of the method, they are carried on MATLAB software with using the convex programming package CVX and the solver SDPT3 4.0. Asteroid 4769 Castalia is used as the target planned for landing, its physical parameters and shape model can be obtained from (Elizabeth and Ludmilla, 2019), the diameter is 1.4km, the density is 2100kg/m3 , the rotation period is 4.095h, and the GM is 94m3 /s2 . The gravity of the asteroid is calculated by the polyhedron gravitation method which is one of the most accurate methods for calculating the gravitational field of asteroids, the shape model used for gravitation calculation includes 2048 vertices and 4092 faces which can be downloaded from (Elizabeth and Ludmilla, 2019). In the simulation, the vehicle is considered to be a rigid body with physical dimensions and its center-of-mass keeps constant despite the fuel consumption. The vehicle parameters are given in Table 1. In the successive solution, the mass mi−1 (k) at each temporal node is calculated by Eq. (41) i−1 with using the thrust Ti−1 (k) and time dilation coefficients λi−1 1 , λ2 , the vehicle’s inertial tensor

Ji−1 (k) used in the ith iteration can be updated by the follow equation with using mi−1 (k). 22

Table 1: Vehicle parameters

Parameters

Values

Initial moment of inertia Mass of vehicle Thrust and Torque Vision sensor parameters Glide-slope Vacuum specific impulse Earths standard gravity constant



  Ji−1 (k) =  

 2940 0 0 2758 0  kg · m2 J= 0 0 0 1974 minit = 1400 kg mdry = 1000 kg Tmax = 20 N Mmax = 0.2 Nm ρS = [0.9, 0, −1.0]T dS = [0, 0, −1.0]T β = 30 deg θ = 15 deg Isp = 225 s g0 = 9.80665 m/s2 

2.10mi−1 (k)

0

0

0

1.97mi−1 (k)

0

0

0

1.41mi−1 (k)

    

(69)

In the simulation, a tolerance of landing error is given, they are 1.5m landing position error, 0.02m/s landing velocity error, 0.01rad/s angular velocity error and 0.02 MRPs error. Moreover, ¯ are needed, a maximum permission iterative number Niter and a iteration cut-off parameter ∆ ¯ the iteration ends and the current solution is the result, here we when i ≥ Niter or k∆i k1 ≤ ∆, ¯ = 10−7 . It must be noted that before starting the successive iteration, a set Niter = 10 and ∆ calculation of the convex subproblem with just considering the dynamics and control constraints is performed, and the obtained solution is used as the initial guess of the reference trajectory. That is to reduce manual selection and ensure that the initial reference trajectory is physically feasible. The initial state and terminal state for the vehicle are given as follows. r0 = [1518.5, −1157.8, −248.93]T m ω 0 = [0, 0, 0]T rad/s

v0 = [−1.423, 1.576, −0.567]T m/s

σ 0 = [0.2897, 0.5303, −0.0895]T

rtf = [682.6374, −132.4631, −89.6853]T m ω tf = [0.001, 0.001, 0.005]T rad/s

vtf = [0, 0, 0]T m/s

σ tf = [−0.1964, −0.3125, 0]T

When the initial reference trajectory and the discrete interval number K1 , K2 remain unchanged,

23

the main factors affecting flight time are the weight of fuel consumption and the size of the ellipsoid enveloping the asteroid. The following will discuss the effect of different  and the three semi-major axes [a, b, c] on the flight time and the trajectory. 5.1. The effect of fuel consumption weight  In this simulation verification, the three semi-major axes for the ellipsoid constraint are set a = 1200 m, b = 750 m, c = 600 m. For discrete interval numbers K1 and K2 , their selection has been pre-estimated and tested, in principle, it is to make the time increment ∆t between the two flight phases not much different. Here, K1 is set to 36 for the descent phase, and K2 is set to 24 for the final landing phase. By definition of the time dilation coefficient, the total flight time is tf = (λ1 + λ2 ) s. Denote ∆m as the fuel consumption, the simulation results of flight time and fuel consumption for different  are given in Table 2 and Fig. 5. Table 2: Flight time and fuel consumption for different 

λ1 λ2 tf (s) ∆m(kg)

=0

 = 10−3

 = 10−2

 = 10−1

301.6580 226.0499 527.7079 12.4905

301.3861 226.5864 527.9725 9.4221

313.3737 233.8214 547.1951 5.8245

410.3230 247.6749 657.9979 3.4450

From table 2 we can know, when  = 0, there is no fuel consumption penalty term in the cost function, the flight time is the smallest, but the fuel consumption is the largest. With the increase of fuel consumption weight, the flight time becomes longer and the fuel consumption decreases, especially when  = 10−1 , the flight time increases dramatically. This indicates that fuel consumption is inversely proportional to flight time over a certain time frame, and the larger , the greater the role of the fuel consumption penalty in the cost function. Therefore, in practical application, the value of  should be reasonably selected according to the thruster’s vacuum specific impulse to ensure that flight time and fuel consumption meet mission requirements. The effect of different  is also reflected in the flight trajectory and the attitude change, the simulation results are shown in Figs. 6-12.

24

Figure 5: Mass m for different 

Figure 6: Trajectories for different 

Figure 7: Displacement r for different 

Figure 8: Velocity v for different 

Figure 9: Thrust T for different 

Figure 10: Angular velocity ω for different 

25

Figure 11: MRPs σ for different 

Figure 12: Torque M for different 

As can be seen from Figs. 7,8,10,11, both the translational motion (r and v) and the rotational motion (ω and σ) change smoothly with the increase of . And from Figs. 9,12, the vehicle requires less fuel and torque to complete the landing as  increases, especially when  = 10−1 , the flight time has increased greatly, the shape of the curve has changed a lot compared to the other three cases. Since the translational motion and the rotational motion are coupled, the thrust direction is controlled by the attitude of the vehicle (as shown in (3), the thrust vector multiplied by the attitude motion related matrix CL S (σ LS )), the less outputs of reaction wheels and the more stable attitude changes, the less transitions on the thrust direction. Therefore, in the application, it is possible to select a large  to meet the flight time and fuel consumption requirements, so that the flight of the vehicle is more stable. 5.2. The effect of ellipsoid size In the simulation, the fuel consumption weight is set  = 10−3 , the discrete interval numbers are still set K1 = 36 and K2 = 24. We only change the magnitude of the three semi-major axes [a, b, c] (m) of the ellipsoid enveloping the asteroid in the ellipsoid constraint, the flight time and fuel consumption are given in Table 3 and Fig. 13. Table 3: Flight time and fuel consumption for different ellipsoid sizes

[a, b, c]

[1000, 650, 550]

[1050, 700, 600]

[1100, 700, 650]

[1200, 750, 650]

λ1 λ2 tf (s) ∆m(kg) Le2g (m)

317.2820 179.7404 497.0224 7.8957 249.3184

315.9930 194.7935 510.7865 8.3160 290.9024

312.2620 206.7027 518.9647 9.1138 323.6111

299.8059 222.2188 522.0247 9.7782 390.7796

26

In Table 3, Le2g represents that distance from the position where the glide-slope constraint begins to the landing site. As the ellipsoid becomes larger, Le2g will gradually increase, from 249.3184m to 390.7796, while the flight time of the descent phase decreases, the flight time of the final landing phase increases, but the total flight time and fuel consumption have increased, especially when a = 1200m, b = 750m, c = 650m, the flight time of the descent phase is less than 300s. The larger Le2g , the smaller flight space of the vehicle, the fewer possible landing paths, the more difficult it is to land, so more fuel is needed. This effect is also reflected in the vehicle’s translational motion and rotational motion, the simulation results are given in Figs. 14-20.

Figure 13: Mass m for different ellipsoid sizes

Figure 14: Trajectories for different ellipsoid sizes

Figure 15: Displacement r for different ellipsoid sizes

Figure 16: Velocity v for different ellipsoid sizes

27

Figure 17: Thrust T for different ellipsoid sizes

Figure 18: Angular velocity ω for different ellipsoid sizes

Figure 19: MRPs σ for different ellipsoid sizes

Figure 20: Torque M for different ellipsoid sizes

In Figs. 13-20, case 1 to case 4 correspond to the 4 different ellipsoid semi-major axes in Table 3, respectively. As can be seen from the Figs. 14-16, both the trajectory curves and the velocity curves have changed, but they have maintained a general trend. In Fig. 17, as the ellipsoid becomes larger, thrusters of the vehicle work at maximum thrust Tmax for more time, so there is more fuel consumption. For the vehicle’s attitude change, the angular velocity and the MRPs curves are almost the same for case 1 to case 3, and the output of reaction wheels is also small. But for case 4, when a = 1200m, b = 750m, c = 650m, ωz and σ3 have produced significant changes after 300s, and reaction wheels work at maximum torque Mmax for most of the time. In summary, increasing the size of the ellipsoid will increase flight time and fuel consumption. Therefore, in practical application, the three semi-major axes of the constrained ellipsoid should be reasonably selected to meet mission requirements.

28

6. Conclusion This paper proposes a trajectory optimization algorithm for 6-DOF asteroid landing trajectory with two consecutive flight phases, since the final time of two flight phases is unknown, we aim to give the optimal flight time and trajectory simultaneously, so this algorithm is called timeoptimal trajectory planning for asteroid landing with two-phase free final time. The algorithm deals with the time-optimal trajectory planning problem under the convex optimization framework. So convexification and discretization are implemented to transform the original continuous-time nonconvex problem into a discrete convex problem. In order to solve the two-phase free final time problem, we normalize the flight time intervals of the two phases into [0, 1] with using two time dilation coefficients λ1 and λ2 , so the free final time problem turns to be a fixed-time problem by minimizing λ1 + λ2 . The cost is also regularized by augmenting a fuel consumption penalty to control the fuel consumption, and the velocity increment is used to give a convex expression of the fuel consumption penalty. Then the successive solution method is used to solve the convex problem and generate the optimal trajectory. Convergence analysis of the method is given, and it is proved that the solutions of a sequence of constrained subproblem will recover the local optimality of the original problem. Meanwhile, the safe-guarding technique of trust regions and virtual control is used to eliminate problems arising from linearization of making the problem unbounded and introducing the artificial infeasibility. Two scenarios for landing on asteriod 4769 Castalia are simulated to analyze the effect of the weight of fuel consumption and the size of the ellipsoid enveloping the asteroid on flight time. The simulation results show that the greater fuel consumption weight, the larger flight time and the less fuel consumption, while the larger ellipsoid, the larger flight time and fuel consumption. So the weight and the ellipsoid size should be reasonably selected according to the vacuum specific impulse of the vehicle and mission requirement. The method can be extended to the free final time problem with more flight phases, and is suitable for application that requires onboard rapid calculation.

Acknowledgments This work was supported by the National Natural Science Fund of China [grant numbers 61503102, 61673057].

29

References Acikmese, B., Ploen, S.R., 2007. Convex Programming Approach to Powered Descent Guidance for Mars Landing. J. Guid. Control Dyn. 30, 1353–1366. Https://doi.org/10.2514/1.27553. Berry, K., Sutter, B., May, A., Williams, K., Barbee, B.W., Beckman, M., Williams, B., 2013. OSIRIS-REx touch-and-go (TAG) mission design and analysis, in: 36th Annual AAS Guidance and Control Conference, AAS, Breckenridge, Colorado. p. 095. Blackmore, L., A¸cıkme¸se, B., Carson III, J.M., 2012. Lossless convexification of control constraints for a class of nonlinear optimal control problems. Syst. Control Lett. 61, 863–870. Https://doi.org/10.1016/j.sysconle.2012.04.010. Blackmore, L., Acikmese, B., Scharf, D.P., 2010. Minimum-Landing-Error Powered-Descent Guidance for Mars Landing Using Convex Optimization. J. Guid. Control Dyn. 33, 1161–1171. Https://doi.org/10.2514/1.47202. Boyd, S., Vandenberghe, L., 2004. Convex optimization. Cambridge university press. Carson, J.M., A¸cikme¸se, B., Blackmore, L., 2011. Lossless convexification of Powered-Descent Guidance with non-convex thrust bound and pointing constraints, in: Proceedings of the 2011 American Control Conference, IEEE. pp. 2651–2656. Crane, L., 2019. The final frontier. Dueri, D., A¸cıkme¸se, B., Scharf, D.P., Harris, M.W., 2016. Customized real-time interior-point methods for onboard powered-descent guidance. J. Guid. Control Dyn. 40, 197–212. Elizabeth,

W.,

Ludmilla,

K.,

2019.

NASA

PDS:

Small

Bodies

Node.

https://pdssbn.astro.umd.edu/. Accessed May 21, 2019. Ferrari, F., Lavagna, M., 2018. Ballistic landing design on binary asteroids: The aim case study. Adv. Space Res. 62, 2245–2260. Https://doi.org/10.1016/j.asr.2017.11.033. Furfaro, R., Cersosimo, D., Wibben, D.R., 2013. tiple sliding surfaces guidance techniques. Https://doi.org/10.2514/1.58246.

30

Asteroid precision landing via mul-

J. Guid. Control Dyn. 36,

1075–1092.

Guo, Y., Hawkins, M., Wie, B., 2013.

Applications of Generalized Zero-Effort-Miss/Zero-

Effort-Velocity Feedback Guidance Algorithm.

J. Guid. Control Dyn. 36, 810–820.

Https://doi.org/10.2514/1.58099. Hu,

H.,

small

Zhu, bodies

S.,

Cui,

with

P.,

reduced

2016.

Desensitized optimal trajectory for landing on

landing

error.

Aerosp.

Sci.

Technol.

48,

178–185.

Https://doi.org/10.1016/j.ast.2015.11.006. Kawaguchi, ence

J.,

Fujiwara,

accomplishment

A.,

Uesugi,

summary

and

T.,

2008.

Hayabusa-2.

Hayabusa-Its technology and sciActa

Astronaut.

62,

639–647.

Https://doi.org/10.1016/j.actaastro.2008.01.028. Lan, Q., Li, S., Yang, J., Guo, L., 2014a. Finite-time control for soft landing on an asteroid based on line-of-sight angle. J. Franklin I. 351, 383–398. Https://doi.org/10.1016/j.jfranklin.2013.08.012. Lan, Q., Li, S., Yang, J., Guo, L., 2014b.

Finite-time soft landing on asteroids us-

ing nonsingular terminal sliding mode control.

T. I. Meas. Control 36,

216–223.

Https://doi.org/10.1177/0142331213495040. Lauretta, D., Balram-Knutson, S., Beshore, E., Boynton, W.V., dAubigny, C.D., DellaGiustina, D., Enos, H., Golish, D., Hergenrother, C., Howell, E., et al., 2017. OSIRIS-REx: Sample Return from Asteroid (101955) Bennu. Space Sci. Rev. 212, 925–984. Https://doi.org/10.1007/s11214017-0405-1. Li, S., Cui, P., Cui, H., 2006. Autonomous navigation and guidance for landing on asteroids. Aerosp. Sci. Technol. 10, 239–247. Https://doi.org/10.1016/j.ast.2005.12.003. Liao-McPherson, D., Dunham, W.D., Kolmanovsky, I., 2016. Model Predictive Control Strategies for Constrained Soft Landing on an Asteroid, in: AIAA/AAS Astrodynamics Specialist Conference, Long Beach, California. p. 5507. Lunghi, P., Lavagna, M., Armellin, R., 2015. A semi-analytical guidance algorithm for autonomous landing. Adv. Space Res. 55, 2719–2738. Https://doi.org/10.1016/j.asr.2015.02.022. Mao, Y., Dueri, D., Szmuk, M., A¸cıkme¸se, B., 2017.

Successive Convexification of Non-

Convex Optimal Control Problems with State Constraints. IFAC-PapersOnLine 50, 4063–4069. Https://doi.org/10.1016/j.ifacol.2017.08.789. 31

Mao, Y., Szmuk, M., A¸cıkme¸se, B., 2016. Successive convexification of non-convex optimal control problems and its convergence properties, in: 2016 IEEE 55th Conference on Decision and Control (CDC), IEEE, Las Vegas, NV, USA. pp. 3636–3641. Mao, Y., Szmuk, M., Xu, X., Acikmese, B., 2018. Successive convexification: A superlinearly convergent algorithm for non-convex optimal control problems. arXiv preprint arXiv:1804.06539 . Pinson, R.M., Lu, P., 2018.

Trajectory Design Employing Convex Optimization for

Landing on Irregularly Shaped Asteroids.

J. Guid. Control Dyn. 41,

1243–1256.

Https://doi.org/10.2514/1.G003045. Reynolds, T., Mesbahi, M., 2017. Small Body Precision Landing via Convex Model Predictive Control, in: AIAA SPACE and Astronautics Forum and Exposition, Orlando, FL. p. 5179. Richards, A., How, J.P., 2002. Aircraft trajectory planning with collision avoidance using mixed integer linear programming, in: Proceedings of the 2002 American Control Conference (IEEE Cat. No. CH37301), IEEE. pp. 1936–1941. Song, Y., Gong, S., 2019. Solar-sail trajectory design for multiple near-earth asteroid exploration based on deep neural networks. Aerosp. Sci. Technol. Https://doi.org/10.1016/j.ast.2019.04.056. Szmuk, M., Acikmese, B., 2018. Successive Convexification for 6-DoF Mars Rocket Powered Landing with Free-Final-Time, in: 2018 AIAA Guidance, Navigation, and Control Conference, Kissimmee, Florida. p. 0617. Szmuk, M., Acikmese, B., Berning, A.W., 2016. Successive Convexification for Fuel-Optimal Powered Landing with Aerodynamic Drag and Non-Convex Constraints, in: AIAA Guidance, Navigation, and Control Conference, San Diego, California, USA. p. 0378. Xiangyu,

H.,

Hutao,

C.,

Pingyuan,

C.,

2004.

tion and guidance for soft landing on asteroids.

An

autonomous

optical

Acta Astronaut. 54,

naviga763–771.

Https://doi.org/10.1016/j.actaastro.2003.09.001. Yang,

H.,

Bai,

X.,

Baoyin,

H.,

2017a.

Finite-time

ing and landing via terminal sliding-mode guidance. Https://doi.org/10.1016/j.actaastro.2016.12.012. 32

control

for

asteroid

hover-

Acta Astronaut. 132, 78–89.

Yang, H., Bai, X., Baoyin, H., 2017b.

Rapid Generation of Time-Optimal Trajectories

for Asteroid Landing via Convex Optimization.

J. Guid. Control Dyn. 40, 628–641.

Https://doi.org/10.2514/1.G002170. Yang,

H.,

Bai,

X.,

Baoyin,

H.,

2017c.

Rapid Trajectory Planning for Asteroid

Landing with Thrust Magnitude Constraint.

J. Guid. Control Dyn. 40, 2713–2720.

Https://doi.org/10.2514/1.G002346. Yang, H., Baoyin, H., 2015. Fuel-optimal control for soft landing on an irregular asteroid. IEEE T. Aero. Elec. Sys. 51, 1688–1697. Https://doi.org/10.1109/TAES.2015.140295. Zexu, Z., Weidong, W., Litao, L., Xiangyu, H., Hutao, C., Shuang, L., Pingyuan, C., 2012. Robust sliding mode guidance and control for soft landing on small bodies. J. Franklin I. 349, 493–509. Https://doi.org/10.1016/j.jfranklin.2011.07.007.

33