Trajectory optimization for lunar soft landing with complex constraints

Accepted Manuscript Trajectory Optimization for Lunar Soft Landing with Complex Constraints Huiping Chu, Lin Ma, Kexin Wang, Zhijiang Shao, Zhengyu Song PII: DOI: Reference:

S0273-1177(17)30534-3 http://dx.doi.org/10.1016/j.asr.2017.07.024 JASR 13334

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

12 May 2016 9 July 2017 14 July 2017

Please cite this article as: Chu, H., Ma, L., Wang, K., Shao, Z., Song, Z., Trajectory Optimization for Lunar Soft Landing with Complex Constraints, Advances in Space Research (2017), doi: http://dx.doi.org/10.1016/j.asr. 2017.07.024

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.

Trajectory Optimization for Lunar Soft Landing with Complex Constraints Huiping Chu1, Lin Ma2, Kexin Wang3 and Zhijiang Shao4,* Zhejiang University, Hangzhou, 310027, People’s Republic of China Zhengyu Song5,* National Key Laboratory of Science and Technology on Aerospace Intelligent Control, Beijing, China, 100854

Abstract A unified trajectory optimization framework with initialization strategies is proposed in this paper for lunar soft landing for various missions with specific requirements. Two main missions of interest are Apollo-like Landing from low lunar orbit and Vertical Takeoff Vertical Landing (a promising mobility method) on the lunar surface. The trajectory optimization is characterized by difficulties arising from discontinuous thrust, multi-phase connections, jump of attitude angle, and obstacles avoidance. Here R-function is applied to deal with the discontinuities of thrust, checkpoint constraints are introduced to connect multiple landing phases, attitude angular rate is designed to get rid of radical changes, and safeguards are imposed to avoid collision with obstacles. The resulting dynamic problems are generally with complex constraints. The unified framework based on Gauss Pseudospectral Method (GPM) and Nonlinear Programming (NLP) solver are designed to solve the problems efficiently. Advanced initialization strategies are developed to enhance both the convergence and computation efficiency. Numerical results demonstrate the adaptability of the framework for various landing missions, and the performance of successful solution of difficult dynamic problems. Keywords: complex constraints; initialization strategies; trajectory optimization for lunar soft landing.

1. Introduction Early in 1959, the lunar exploration began when the Soviet Union proposed the Luna program. The Luna 2 was launched. This was the first spacecraft to reach the lunar surface and successfully made the first crash landing on the Moon. The United States began the Apollo program in 1961, which accomplished the transporting of humans to the Moon [1]. Recently, studies for lunar exploration have been pursued in many countries [2]. Chinese lunar exploration program, which began in 2004, is divided into three main operational phases, including orbital missions, soft landing mission, and sample return mission. The Chandrayaan-2 developed by the Indian Space Research Organization has a planned launch to the Moon by the end of 2017 or by the beginning of 2018. This will be India's second lunar exploration mission after Chandrayaan-1. Luna-Glob, a lunar exploration program by the Russian Federal Space Agency, will be the first of the missions planned before the creation of a fully robotic lunar base, and has a planned launch on 2024. The Japanese Lunar Exploration Program is a project of robotic and human missions to the Moon, and its main goal is to elucidate the origin and evolution of the Moon and utilize the Moon in the future. After the Apollo * Correspongding author. Email addresses: [email protected] (H. Chu), [email protected] (L. Ma), [email protected] (K. Wang), [email protected] (Z. Shao), [email protected] (Z. Yu) 1 Graduate Student, College of Control Science and Engineering, Zhejiang University. 2 Ph.D. Candidate, College of Control Science and Engineering, Zhejiang University. 3 Research Scientist, College of Control Science and Engineering, Zhejiang University. 4 Professor, College of Control Science and Engineering, Zhejiang University, State Key Laboratory of Industrial Control Technology. 5 Vice Chief Designer, National Key Laboratory of Science and Technology on Aerospace Intelligent Control.

Program from 1961 to 1972, NASA officials have not lost interest in the Moon. NASA has recently launched Lunar CATALYST, which is a program designed to assist and encourage private companies that are interested in lunar exploration [3]. Soft landing is of great importance to lunar exploration, and is also one of the primary techniques for lunar exploration programs of all the abovementioned countries. For different types of lunar soft landing missions, safety, precision, cost, and efficiency are the key issues. The studies on trajectory optimization for lunar soft landing can effectively contribute to improving these issues. Especially, China planned to land softly on the back side of the Moon in 2018, which will become the first country to get there. The back side of the Moon is the hemisphere of the Moon that always faces away from Earth, and it is always shielded from radio transmissions from the Earth. In addition, the terrain of far side is rugged. Therefore, improving the autonomy and intelligence of the lander becomes a real demand. Two typical main missions of the lunar soft landing are studied in this paper, including Apollo-like Landing mission from low lunar orbit and Vertical Takeoff Vertical Landing (VTVL) mission from lunar surface. The Apollo-like Landing mission is a typical descent and landing process from the lunar orbit by the Apollo descent planning [4]-[5][35], as shown in Fig. 1. After the lander and the Command and Service Module (CSM) have separated to a safe distance in the Parking Orbit (approximately 100 km altitude), the lander performs the Descent Orbit Initiation (DOI) and is inserted to the Hohmann transfer orbit. The purpose of the Hohman Transfer is to efficiently reduce the altitude with no power consumed from approximately 100 km to 15 km (the perilune altitude) in preparation for the Powered Descent Initiation (PDI). Thereafter, the basic trajectory design for the powered descent is established as three phases to satisfy the operational requirements. These three phases are called the braking phase, the approach phase, and the landing phase. The braking phase follows the PDI starting at perilune altitude, and is designed primarily for the orbit velocity being reduced, during which the altitude is decreased from approximately 15 km to 2.2 km. The approach phase is designed for the pilot with visibility of the landing area and for monitoring the approach to the lunar surface. The altitude is reduced from approximately 2.2 km to 152 m. Finally, the basic purpose of the landing phase is to provide a flight at low velocity and pitch attitude perpendicular to the lunar surface. The attitude angle should be a smooth transition among the three phases. The initial trajectory analysis which led to this design was performed by Bennett and Price [4]. A detailed description of the lunar lander landing design strategy was provided by Cheatham and Bennett [5]. Powered lunar landing is studied in this paper. It was reported that the descent engine nominally provides the upper fixed position of the throttle about 43148 N of thrust (maximum thrust) and is always throttleable in the region of approximately 4671 N of thrust to 28023 N (65% of maximum thrust) of thrust [6]-[7]. A gap in the data appears because the requirement for throttling ranges from 28023 N of thrust to maximum thrust as a result of a serious erosion problem in that range [6]. Thus, the thrust magnitude cannot be adjusted continuously between maximum and minimum values, thereby leading to difficulties in solving the optimal trajectory.

Fig. 1 Lunar Landing Phase for Apollo Descent Strategy

Another mission is called the Vertical Takeoff Vertical Landing (VTVL) [8], which might be a promising lander mobility method both for the Earth and the Moon. As of 2016, VTVL is under intense development as a promising technology for reusable rockets, with two companies, Blue Origin and SpaceX, both having demonstrated recovery of launch vehicles, with Blue Origin's New Shepard booster rocket making the first successful vertical landing following

a test flight that reached outer space, and SpaceX's Falcon 9 Flight 20 marking the first landing of a commercial orbital booster. Since the Moon is surrounded by vacuum, this might be the only mobility method for long distance travel on the surface for lunar exploration in different regions. Such a VTVL lander would take off vertically, translate a desired distance, and vertically land softly onto the lunar surface, as shown in Fig. 2. Several reasons are listed for using this mobility method with the use of lunar rovers and robots. First, this mobility method can handily explore more complex areas on the lunar surface because of the ability to cross different terrains, such as craters and boulder fields. Second, with the creation of lunar base in the future, the mobility method of VTVL contributes to lunar surface activities among lunar bases and exploration destinations. Because of the constraints on landing attitude angle, attitude angular rate, and velocity, the vertical landing is significantly more difficult compared with the vertical takeoff. In a previous study [8], short distance, low height trajectory of VTVL was optimized with two-dimensional dynamics model. The dynamics model is extended in this paper with three-dimensional coordinates that remain unified with the Apollo-like Landing mission. Approximately 83% of the Moon’s surface is composed of lunar highlands [9]. This study aims at the trajectory optimization design of avoiding terrain obstacles during Vertical Takeoff Vertical Landing in particular.

Fig. 2 The Flight of Vertical Takeoff Vertical Landing (VTVL)

The trajectory optimization problem is generally converted to the Optimal Control Problem (OCP). Two general categories for solving OCP exist, namely, indirect and direct techniques [10]. The indirect method based on the Pontryagin Maximum Principle converts the optimal control problem to a Hamilton boundary-value problem (HBVP) [11]-[12]. The solution solved by indirect method with high accuracy satisfies the first-order optimality conditions. However, several disadvantages exist, including small radii of convergence, the potential difficulty in formulating the HBVP, an accurate initial guess for the costate, and knowledge of the switching structure for path-constrained problems. In a direct method, the optimal control problem is discretized and transcribed directly into a nonlinear programming problem (NLP) without formulating an alternate set of optimality conditions [13]. Compared with the indirect method, the direct methods have advantages as follows: deriving the first-order necessary conditions is not necessary; no demand for an accurate initial guess; and no requirement for an initial guess for the costate. This paper is based on the Gauss Pseudospectral Method (GPM), which has advantage over other numerical methods because of various characteristics, including high accuracy in both the primal and dual solutions and equivalence between the direct and indirect forms [14]-[20]. The optimization problem of spacecraft trajectories based on GPM has been studied extensively [21]-[24]. The optimal reconfiguration of spacecraft formations was studied in [21], the results of which demonstrate the applicability of the GPM to spacecraft trajectory optimization problems. Constrained trajectory optimization using GPM is presented in [22], thereby demonstrating that GPM is capable of providing highly accurate solutions. The research in [23] is based on the successful determination of GPM and verification of a numerical technique that can generate an optimal reentry trajectory that satisfies all of the specified vehicle and waypoint and no-fly zone constraints. In [24], the offline reentry trajectory is designed based on GPM, which is used as the initial value guess for real-time reentry trajectory optimization. The lunar soft landing has been investigated for many years, but most papers focus on the specific landing scenarios [25]-[27]. However, the requirements and constraints cause the diversity of lunar soft landing missions. These diverse landing missions are expected to be resolved in a unified trajectory optimization framework. The two abovementioned typical main missions of the lunar soft landing are investigated with unified optimization method in this paper. Their three-dimensional dynamics model and cost function are presented in a unified way, and their constraints are described

with similar structure. In particular, the discontinuous adjustable thrust magnitude constraint is described as a path constraint by using R-function method [28]-[29].The unified optimization framework proposed in this paper is robust and adaptable for solving diverse landing missions, and flexible for handling complex constraints. However, the initial values and solution convergence problems should be concerned. Several initial guess strategies in this paper can effectively enhance solution convergence and contribute to improving accuracy in obtaining an optimal trajectory. 2. Lunar lander formulation Two main missions of the lunar soft landing are discussed in this paper, including the Apollo-like Landing and the VTVL. Even though the two missions differ in requirements and constraints, their three-dimensional dynamics model and cost function are the same. It is expected that the different lunar soft landing missions could be optimized in a unified way. In this section, we will first discuss the model assumptions. Then, we will provide the coordinate systems and unified dynamic model, as well as objective function. 2.1. Assumptions To capture the limitations of real spacecraft while maintaining an adaptable model, several simplifications are made. The Moon is assumed to be spherical, and its rotation is ignored. Besides, the lander is modeled as a point mass in a uniform gravitational field with the thrust acting directly on the center of mass. While a point mass does not have a meaningful attitude, the lander’s body axes is assumed to coincide with the thrust vector. Therefore, the attitude angles of the lander are defined by the thrust direction angles, which are called pitch and yaw angles. 2.2. Coordinate Systems and Lander Model In [30], the reentry spacecraft dynamics equation was derived, using three variables to describe speed parameters, including velocity magnitude, flight path angle, and azimuth. However, the reciprocal of the velocity, 1/V , is in the motion equations. When the velocity of lander is reduced to zero during lunar soft landing trajectory optimization, the solution will be singular. To overcome this difficulty, a reasonable three-dimensional dynamics model is given. Two right-handed reference systems (illustrated in Fig. 3) are defined, as follows: (1) OXYZ is the lunar central inertial coordinate, where Y-axis points to the longitude of 0◦ in the lunar equatorial plane, and X-axis points to the arctic pole of the lunar. (2) Axyz is the motion coordinate, where the x-axis is along to the direction of OA, and z-axis is located in the AOX plane and perpendicular to Ax.

Fig. 3 Coordinates of Lunar Soft Landing

In Fig. 3, r =(r, α, β) is the position vector. r is the radial distance from the lunar center to the lander. The longitude α is positive to the East. The latitude β is positive to the North, and the latitude of 0◦ is located in the lunar equatorial plane. A is the position of the lander. T =(T, ϕ, ψ) is the thrust vector. T is the thrust magnitude. ϕ and ψ are the pitch and yaw angles, respectively. Based on the Newton's second law, dynamic equations are derived and given as follows:





r  Vx，   

Vx  

Vy  

Vz 

T sin 



m T cos  sin 

Vz r

V V 2 y

T cos  cos

m

，

r cos 

g

m



Vy



2 z

r VxVy



r

VxVz r



VyVz tan 

(1)

r

V tan  2 y

r

where Vx, Vy , and Vz are the velocities in the motion coordinate. m is the mass of the lander. g=-μ/r2 is the Moon gravitational acceleration, where μ is the Moon gravitational coefficient. The mass of the lander varies with time and equation for the mass flow rate is as follows: 

m   T I sp g

(2)

where the effective specific impulse of the lander, Isp, is assumed to be a constant across the entire throttle range, and g⊕ is the Earth gravitational acceleration at sea level. Additionally, the allowable attitude angular rate should be constrained to avoid sudden change of the attitude angle, which is expected to shift smoothly during the flight. Therefore, two state equations are introduced accordingly as follows: 

  ，



  

(3)

where ωϕ and ωψ are the attitude angular rate, called pitch rate and yaw rate, respectively. The purpose of the optimal control strategy is to guide the lander from the initial to the final conditions with the minimum amount of fuel consumption, thereby leading to the so-called minimum-fuel lunar soft landing problem with the following cost function.

J  m  t f 

(4)

3. Design of Missions for Lunar Soft Landing The three-dimensional dynamics model and cost function are the same in two typical landing missions (Apollo-like Landing and VTVL). According to the specific requirements and constraints, these landing missions are both designed to several scenarios to adapt to different situations and applications. 3.1. Apollo-like Landing Two scenarios are designed for the Apollo-like Landing as follows: the Three-Phase Landing from Pinpoint (TPLP), and the One-Phase Landing from the Orbit (OPLO). The requirements and constraints of the TPLP mission are same as the Apollo 12 landing plan. The entire process of soft landing is divided into three phases, including the braking phase, approach phase, and landing phase, as shown in Fig. 1. The altitude of the braking phase decreases from approximately 15 km to 2.2 km. In addition, at the final time of the braking phase, the horizontal distance to the landing target position on the lunar surface, as well as the vertical and horizontal velocities, should be lower than specific values. The approach phase reduces the height from 2.2 km to 152 m, and vertical and horizontal velocities at the final time are also lower than specific values. Moreover, the attitude is approximately vertical to the lunar surface. The lander is required for soft landing on the surface at the final time of the landing phase, in which velocity is reduced to zero, and attitude is vertical to the lunar surface. For the OPLO mission, the requirements and constraints of the braking phase and the approach phase are eliminated. From the Hohmann transfer orbit rather than from a pinpoint, it lands directly to target landing position on the lunar surface, with velocity reduced to zero and attitude vertical to the lunar surface at final time. Simultaneous

consideration of the trajectory of the lander and its initial descent position on the orbit gives more freedom to optimizing the landing mission. Here the transfer orbit is specified by five orbital elements, including semi-major axis, eccentricity, orbital inclination, longitude of the ascending node, and argument of perilune. It can be expected that OPLO is more fuel-efficient than TPLP. 3.2. Vertical Takeoff Vertical Landing VTVL is a mobility method on the lunar surface. Such a VTVL lander would take off vertically, transfer a desired horizontal distance, and land softly on lunar surface. The attitudes at both the initial and final time are constrained to be vertical to the surface. There is one scenario for the VTVL studied in this paper, namely the constrained VTVL with obstacle avoidance which means the obstacles must be considered explicitly during the flight.

4. Constraints for Lunar Soft Landing The two lunar landing missions have constraints with similar structure, although they still have their specific ones. All of these constraints are described below, including the boundary constraints, the inequality and equality path constraints, and the checkpoint constraints. 4.1. Boundary Constraints Boundary constraints include the initial and terminal conditions. Generally, for the TPLP and VTVL missions, the initial conditions refer to the initial position in Eq. (5), the initial velocity in Eq. (6)，the initial mass in Eq. (7), and the initial attitude in Eq. (8), at initial time t0, as follows: (5) r  t0  =Rmoon  h0 ,   t0  = 0 ,   t0  =0

Vx  t0  =Vx ,0 , Vy  t0  =Vy ,0 , Vz  t0  =Vz ,0

(6)

m  t0  =m0

(7)

  t0  =0

(8)

where RMoon is the Moon average radius. h=r-RMoon is the altitude to the lunar surface. And the subscript ‘0’ denote the initial state at t0. For elliptical orbit, semi-major axis (a) and eccentricity (e) describe the form of the orbit. Orbital inclination (i), longitude of the ascending node (Ω), and argument of perilune (ω) determine the direction of the orbit in outer space. The true anomaly (υ) defines the position along the orbit. The six elements are always uniquely calculated by the position and velocity vector, which are converted by orbital mechanics formula [31]. Five elements, including a, e, i, Ω, and ω are constant and determine the initial elliptical orbit. The true anomaly will be optimized, which will determine the initial descent position of lander on the orbit. For the OPLO mission, the initial conditions refer to the orbital elements in Eq. (9), the initial mass in Eq. (7), and the initial attitude in Eq. (8).

a  t0  =aset ,

e  t0  =eset，

  t0  =set ,

  t0  =set

i  t0  =iset ,

(9)

The terminal condition refers to the position and state of the lander at final time tf. They may be equality or inequality constraints. The subscript ‘f’ denote the final state at tf. The terminal altitude is given as follows:

r  t f  =Rmoon  h f

(10)

The horizontal distance to the landing target position on the lunar surface at final time is given as follows:

s f ,min  s  t f   s f ,max

(11)

where

s  t f  =    f  Rmoon       f  Rmoon  2

2

And S is the distance in the horizontal plane between two locations described by longitude and latitude. The terminal longitude and latitude are as follows:

  t f  = f ,   t f  = f

(12)

The terminal velocities are given as Eq. (13)or Eq. (14)

V, f  Vx  t f  =Vx , f V||, f 



Vy  t f

   V t   2

z

(13)

2

f

V, f ,min  V, f  Vx  t f   V, f ,max V||, f ,min  V||, f 

V t   V t  2

y

f

z

2

f

 V||, f ,max

(14)

where the V and V|| refer to the vertical velocity and horizontal velocity. The terminal attitude is constrained as Eq. (15) or Eq. (16)

  t f  = f

(15)

 f ,min    t f    f ,max

(16)

4.2. Inequality and Equality Path Constraints For the path constraints, the variables should satisfy the upper or lower bound values or specific requirements during the flight. All the subscripts ‘max/min’ denote the maximum/minimum allowable values throughout the paper. First, the lander should fly over the lunar surface and under the specific altitude.

0  h  hmax  Rmoon  r  Rmoon  hmax 

(17)

Second, the attitude angle constraints are given as follows:

min    max ,  min     max

(18)

Third, to have smooth angle change, attitude angular rates are constrained as follows: 



 =    ,max ,  =    ,max

(19)

Fourth, the descent engine nominally provide upper fixed position of the throttle (maximum thrust) and is always throttleable in the region of approximately minimum thrust to 65% of maximum thrust [6]-[7]. This means that the thrust magnitude cannot be adjusted continuously in full range between upper and lower values. It is formulated as follows: Tmin  T  65%*Tmax or T  Tmax (20) The discontinuity of the thrust magnitude constraints with two separated ranges is the main difficulty for optimizing the trajectory. With introduction of R-function, the Boolean operation of set can be converted to an algebra operation of function [28]-[29]. Union: f1

f 2  f1  f 2 + f12  f 2 2

(21)

For thrust magnitude constraints, we can have the following definitions:

f1  (T  Tmin )*(65%* Tmax  T ) f 2  (T  Tmax )*(Tmax  T )

(22)

and the following R-function:

F  f1

f 2  f1  f 2 + f12  f 2 2  0

(23)

where F is the union of the separate sets. With path constraint F  0 , the thrust magnitude is constrained in the range from Tmin to 65%Tmax, or just fixed at the upper value.

In addition, for a variety of missions of VTVL, it may be desirable to avoid obstacles like Lunar mountains or boulder fields. The mathematical model of the topography is

      R  2       R  2  0i moon  0i moon    Rmoon  h h  ,     hi  exp        ti ti i   M

(24)

Where α and β defer the longitude and latitude on the horizontal plane of the Moon, and h is the height variable of terrain. ∆h is used to let the minimum height of lunar terrain be less than zero. hi ，αti ，βti are the parameters that determine the terrain profile. i represents the i-th mountain, and M is the total number of Lunar mountains. The simulation of Lunar topography is shown in Fig. 4. For safety consideration, margins are specified as ∆hmargin.

h  h  ,    hmargin

(25)

Fig. 4 Simulation of lunar terrain

4.3. Checkpoint Constraints The TPLP mission is established as a three-phase maneuver to satisfy the operational requirements. The three phases are linked by two checkpoints. State variables cannot change suddenly at checkpoints. Therefore, the terminal state of the previous stage is the initial state of the next stage. At checkpoint, continuity constraints must be satisfied as follows:

rt   =r  p

p 1

f

 p

Vx ,t =Vx

 t0  ,  t p  =  p 1  t0  , t p  =  p 1  t0  f

 p 1

f

mt  =m p

f

p 1

f

 t 0  , Vy , t

 p f

=Vy

 p 1

 t0  , Vz,tp  =Vz  p 1  t0  f

(26)

 t0  , t p  =   p 1  t0  , t p  =  p 1  t0  f

f

where the superscript ‘p’ is the p-th phase of lunar soft landing. 5. Optimal Control Problem for Lunar Soft Landing The trajectory optimization problem is an optimal control problem. With the introduction of formulation in Sec. 2, the design of landing missions in Sec.3, and the description of constraints in Sec.4, three soft landing sub-missions described as optimal control problems are presented in Table 1. It is expected that the three scenarios listed in Table 1 are handled within the same trajectory optimization framework. Generally, the following optimal control problem in Mayer form [12] on the time interval t  [t0, tf] is considered. We determine the state x (t), control u (t), initial time t0 , and final time tf , which minimize the following cost

function,



J   x  t0  , t0 , x  t f  , t f

min



(27)

subject to the dynamic equation constraints, the boundary constraints, and the inequality and equality path constraints, as follows: 

x  f  x  t  , u  t  , t  t  t0 , t f 



(28)

 x  t0  , t0 , x  t f  , t f =0



(29)

C  x t  , u t  , t   0

(30)



The optimal control problem is to find optimal control law u (t), which minimizes the performance measure of Eq. (27) while satisfying Eqs. (28)–(30). The optimal control problem of Eqs. (27)–(30) is referred to as the continuous Mayer problem. Table 1 Summary of Three Soft Landing Sub-missions for Optimal Control Problem Differences with other missions Missions

Cost function

Dynamic equation

Braking phase

Apollo-like Landing

Three-Phase Landing from Pinpoint (TPLP)

Approach phase

Eq. (4)

Eq. (1) Eq. (2) Eq. (3)

One-Phase Landing from the Orbit (OPLO)

Eq. (4)

Eq. (1) Eq. (2) Eq. (3)

Constrained VTVL with obstacle avoidance

Eq. (4)

Eq. (1) Eq. (2) Eq. (3)

Landing phase

VTVL

Constraints Boundary constraints Path constraints Initial Terminal conditions conditions Eq. (5) Eq. (10) Eq. (6) Eq. (11) Eq. (7) Eq. (14) Eq. (10) Eq. (17) Eq. (14) Eq. (18) Eq. (16) Eq. (19) Eq. (10) Eq. (23) Eq. (12) Eq. (13) Eq. (15) Eq. (10) Eq. (17) Eq. (12) Eq. (18) Eq. (7) Eq. (13) Eq. (19) Eq. (9) Eq. (15) Eq. (23) Eq. (17) Eq. (5) Eq. (10) Eq. (18) Eq. (6) Eq. (12) Eq. (19) Eq. (7) Eq. (13) Eq. (23) Eq. (8) Eq. (15) Eq. (25)

Checkpoint constraints

Eq. (26)

6. Optimization Method Gauss Pseudospectral Method (GPM) with advanced initialization strategies is utilized to solve the proposed trajectory optimization for lunar soft landing. It is expected that the three scenarios listed in Table 1 are handled within the unified trajectory optimization framework. 6.1. Gauss Pseudospectral Method GPM is a direct transcription method that discretizes and transcribes the continuous Mayer optimal control problem by Eqs. (27)-(30) into a NLP problem, which can then be solved by well-developed NLP solvers such as SNOPT [32]. GPM approximates to the state and control profiles with interpolation polynomials on Gauss collocation points. Advantages of GPM in solving optimal control problems include that the KKT conditions of the NLP are exactly equivalent to the discretized first-order necessary conditions of the optimal control problem [14]-[16]. Moreover, GPM is applicable to a large variety of problems including those with path constraints, boundary constraints and free time. At last, GPM can find accurate approximate solution to the optimal control problem and converge quickly. 6.2. Unified Trajectory Optimization Framework with Initialization Strategies The three lunar landing scenarios described in Sec. 5 are to be handled in the unified optimization framework shown

in Fig. 5. The underlying optimal control problems in the form of Eqs. (27)-(30) are discretized and solved by GPM. And the NLP solver is SNOPT, which is an implementation of a Sequential Quadratic Programming (SQP) algorithm for general nonlinear optimization. Finally, the optimized control and state profiles are solved by SNOPT. However, this framework usually requires solution of a large-scale and nonconvex NLP problem. Special optimization strategies are necessary for reliable and fast convergence, and therefore are crucial to successful implementation of the framework. In particular, the discontinuity of the thrust magnitude constraints Eq. (20) is a challenge to efficient numerical optimization methods as they are commonly designed based on derivative information. Even if the constraints are reformulated by Eqs. (22)-(23), the resulting problem is not easy to converge because of the nonconvexity. A heuristic is applied to treat such thrust domain as continuous for a while, imposing only lower and upper bounds to cover the whole domain to replace the original discontinuous constraints. The bounds are much easier to be handled in optimization and often lead to successful solution. The solution can be infeasible to the original domain, while it usually provides a good starting point to improve convergence of the problem subjected to the original more complex constraints. For the VTVL mission, an obvious and useful strategy is to take the corresponding obstacle-free trajectory as the initial point when dealing with case requiring obstacle avoidance. The size of the NLP problem increases with the number of collocation points in discretization. Finer discretization can approximate the dynamic problem better, while generating a larger NLP problem challenging the performance of NLP solvers. On the other hand, if we reduce the size of the NLP problem by coarse discretization, we may end up with a solution of dissatisfactory accuracy. In order to obtain an accurate approximation to the solution of the dynamic problem efficiently, a NLP problem based on coarse discretization (i.e. 10 Gauss points) is solved firstly. Not much initial values need to be provided for the optimization at this stage, and the convergence is not sensitive to the initial values either. Taking advantage of this solution, for example, by linear interpolation, we can provide initial guess for variables of a larger NLP problem deriving from refined discretization. The initialization can be expected to capture certain characteristics of the profiles we aim to benefit the subsequent solution. Moreover, the total computation time by this refinement strategy can be much less than that of solving a large problem from a cold starting point. Unified trajectory optimization framework for lunar soft landing

Normalize a trajectory optimization problem for lunar soft landing in the form of Eqs.(28)-(31)

Initial-guess-based Iterative solution

Discretization (GPM)

Constraints enhancement

A finite-dimensional NLP problem

Collocation points refinement Advanced Initial guess

Nonlinear Optimization （SNOPT）

Optimized Control profile Optimized State profile Fig. 5 Unified Trajectory Optimization Framework

For the two troubles mentioned above, the advanced initial strategies are called constraints enhancement and collocation points refinement respectively. Ingenious strategies to gain the initial guess are expressed in the following calculation steps. Step 1: Simplify the optimal control problem. Two measures are applied to simplify the optimal control problem. The first effective measure is to take the value of thrust magnitude in the full range of upper and lower bounds; in this manner, the optimal control problem is simplified with simpler path constraints. Another effective measure is to eliminate the troublesome constraints on the base of original optimal control problem. For the constrained VTVL mission with an obstacle, the constraints of the obstacles are not introduced. Furthermore, all of the constraints that

cause difficulties can be considered for elimination to simplify the original optimal control problem. Step 2: Solve the simplified optimal control problem with small collocation points. With small collocation points, such as 10 collocation points, less initial guess values of the variables are required. In this case, the convergence is not sensitive to the initial values. The approximate feasible solution can be calculated by using the unified trajectory optimization framework. Step 3: Solve the simplified optimal control again with more collocation points based on the results in Step 2 as the initial guess values. The control and state profiles calculated in Step 2 are interpolated to the same number of collocation points and used as the initial guess. Similarly, the approximate optimal solutions are obtained on the basis of the unified trajectory optimization framework. Step 4: Solve the original optimal control problem with the same collocation points in Step 3 by using the calculation results in Step 3 as the initial guess values. The control and state profiles calculated in Step 3 are directly used to the initial guess. The accurate optimal trajectory can be gained by solving the original optimal control problem with the unified trajectory optimization framework. In summary, the unified trajectory optimization framework has a total of three optimization calculations in Steps 2, 3 and 4, respectively. And the two calculations in Step 2 and Step 3 are to provide an excellent initial guess values. The excellent initial guess used in the original optimal control from Step 3 can effectively enhance solution convergence, help improve the accuracy for optimal trajectory, and spend less computation time for large-scale problems. 7. Numerical Results In this Section, results are presented to validate the proposed trajectory optimization framework for lunar soft landing, as well as the initial guess strategies. Three scenarios for optimal control problem described in Sec. 5 are optimized by the proposed framework. State variables in the formulation of the dynamic problems include r, α, β, Vx, Vy, Vz, m, ϕ, and φ. Control variables include T, ωϕ, ωψ. The basic parameters are listed in Table 2, and the bounds on variables are given in Table 3. Note that thrust T between its bounds is discontinuous and is reformulated by Eqs. (22) -(23) for optimization. GPM implements discretization based on 45 Gauss collocation points for all the numerical simulations, unless a description explains that a specific discretization strategy is used. Table 2 Basic Parameters for Trajectory Optimization Symbol RMoon g⊕

Value 1738km 9.8m/s2

Symbol μ Isp

Value 4.902802627×10 12m3/s2 302.39s

Table 3 Bounds on Variables Variable h ϕ ψ

Min 0 -20° -180°

Max 30km 90° 180°

Variable ωϕ ωψ T

Min -8°/s -8°/s 4671N

Max 8°/s 8°/s 43148N

7.1. Three-Phase Landing from Pinpoint (TPLP Mission) The initial conditions of the braking phase and the terminal conditions of each phase are listed in Table 4 and Table 5, respectively. The data therein come from the Apollo Landing reports [1][33]-[35]. Moreover, simultaneous optimization of all the three phases needs the checkpoint constraints in Eq. (26) to guarantee phase connection. The numbers of Gauss collocation points for phases braking, approach and landing are 45, 15 and 15, respectively. Table 4 Initial Conditions (TPLP Mission) Variable h0 Vx,0 m0

Value 15700m -7m/s 15234kg

Variable α0 Vy,0

Value -1.43° -1651.9m/s

Variable β0 Vz,0

Value -8.43° 375.3m/s

Table 5 Terminal Conditions (TPLP Mission) Phase Braking phase

Variable h(1) f s(1) f

Value 2.2km ≤9km

Variable (1) |V⊥,f | |V(1)‖,f |

Value ≤50m/s ≤300m/s

h(2) f (2) |V⊥,f | h(3) f α(3) f β(3) f

Approach phase Landing phase

≤4.5m/s ≥70° 0m/s 0m/s 90°

|V(2)‖,f | ϕf(2) (3) V⊥,f V(3)‖,f ϕ(3) f

152m ≤16m/s 0m -23.45° -2.94°

The TPLP mission is calculated in the unified trajectory optimization framework. And the initial guess values are generated by the constraints enhancement and collocation points refinement strategies. The iterative calculation procedures are shown as follows: Step 1: Solve the simplified optimization problem with small collocation points. Based on the original optimization problem, the thrust constraints are imposed only lower and upper bounds to cover the whole domain to replace the original discontinuous constraints. The numbers of Gauss collocation points for three phases are 15, 5 and 5, respectively. The initial guess values for the simplified optimization problem are obtained by the measure of linear interpolation in Table 6, and shown in Fig. 6 with pink circle curves. The coarse optimized trajectories are illustrated in Fig. 6 with red circle curves. The calculation time is 5.58 seconds. Step 2: Solve the simplified optimization problem again with finer collocation points. The coarse optimized trajectories obtained in Step 1 are interpolated to 45 collocation points as the initial guess values. The trajectories for the simplified problem are presented in Fig. 6 with yellow curves, and the optimization solution take 29.58 seconds. Step 3: Solve the original optimization problem with finer collocation points. The original optimization problem is solved with the results in Step 2 as the initial guess values, and the optimized trajectories are illustrated in Fig. 6, where the black curves, green curves, and blue curves refer to the braking phase, approach phase, landing phase respectively. Besides, the calculation time is 25.48 seconds. Table 6 Initial guess values by linear interpolation Variable

Values

r0 

r

β

0 

Vy

Vy 0 

m ψ

r

0

 rf  t

Variable α

0 

Vx

Vx 0 

Vz

Vz 0 

tf

 V

0

  f t tf

y0

 Vy f  t tf

ϕ tf

m0 0°

Values



0

 f t tf

V

x0

 Vxf  t tf

V

z0

 Vz f  t tf

0° 800s

0 -5

Longitude ,deg

Altitude H,m

15000

10000

5000

-10 -15 -20 -25

0 0

200

400 Simulation Time,s (a)Altitude H

600

800

-30

0

200

400 Simulation Time,s (b)Longitude 

600

800

0 0

Velocity Vx,m/s

Latitude ,deg

-2 -4 -6

-40

-60

-8

-80

-10

0

200

400 Simulation Time,s (c)Latitude 

600

800

0

0

200

400 Simulation Time,s (d)Velocity Vx

600

800

0

200

400 Simulation Time,s (f)Velocity Vz

600

800

0

200

400 Simulation Time,s (h)Yaw 

600

800

400

-500

Velocity Vz,m/s

Velocity Vy,m/s

-20

-1000

300 200 100

-1500 0 0

200

400 600 Simulation Time,s (e)Velocity Vy

800

100

0

80

-5

60

-10

Yaw  ,deg

Pitch ,deg

-2000

40 20 0 -20

-15 -20 -25

0

200

400 Simulation Time,s (g)Pitch 

600

800

-30

4

5

x 10

Thrust T,N

4 3 2 1 0 0

200

400 Simulation Time,s (i)Thrust T

600

800

Fig. 6 Optimized Landing Trajectories (TPLP Mission)

If the initial guess strategy of constraints enhancement is not used in solving the optimization problem, it fails to obtain solution of the TPLP mission. If the collocation points refinement strategy is not adopted, it takes 206.62 seconds to gain the optimized trajectories of simplified problem. However, the iterative calculation procedures in Step 1 and Step 2 take only 35.16 seconds in total due to the application of the collocation points refinement strategy. Therefore, the two ingenious strategies to generate the initial guess values can effectively enhance solution

convergence, and sufficiently improve computation efficiency. The duration of this three phases landing problem is tf =637.15s. And the fuel consumption is 6962.43kg. Fig. 6 (a)-(c) present the positon of lander vs. time. The Fig. 6 (a)-(c) show that the height of lander decreases gradually and the lander arrives two checkpoints. Finally, the lander land in the desired position at (0m, -23.45°, -2.94°). Fig. 6 (d)-(f) shows that the velocities are reduced to zero eventually. In addition, the attitude angles in Fig. 6 (g)–(h) are smooth during the flight, as the result of imposing attitude angular rate constraints Eq. (19). Therefore, radical changes in ° attitude angles due to the final condition ϕ(3) f =90 are avoided. Fig. 6 (i) gives the optimal thrust, which satisfies the discontinuous magnitude constraints Eq. (20). Along the optimal paths, the thrust always switches from the two extremes (minimal or maximal) values. 7.2. One-Phase Landing from the Orbit (OPLO Mission) The initial conditions are given by Table 7, where the orbital elements follow Apollo 12 landing reports [34]. And the terminal conditions are the same to the Landing phase of Sec. 7.1 in Table 5. The checkpoint constraints are not applicable to this mission. Table 7 Initial Conditions (OPLO Mission) Variable aset Ωset

Value 1801369.46m 325.7370°

Variable eset ωset

Value 0.0269° 335.3317°

Variable iset m0

Value 164.7122° 15234kg

The unified trajectory optimization framework as well as the application of initial guessing strategies are used to calculate the OPLP mission. The calculation procedures are the same as that of Sec. 7.1. The numbers of collocation points for three phases are 15 firstly. The initial guess values of the simplified optimization problem and the coarse optimized trajectories are shown in Fig. 7 with pink circle lines and red circle curves respectively, and 1.92 seconds were taken during this simulation. Making the interpolated red curve as the initial guess values, the optimization solution of fine optimized trajectory for simplified problem takes 6.65 seconds. Since the thrust profiles does not violate the discontinuous thrust constraints, no further solution is needed. Profiles of the states and controls are shown by blue curves in Fig. 7. The thrust profile in Fig. 7 (i) satisfies the discontinuous magnitude constraints Eq. (20). Along the optimal paths, the thrust has been at the maximum in most cases and was throttleable between minimal thrust and 65% of full thrust at the end of the flight. In addition, the calculation time without application of collocation points refinement strategy is 18.62 seconds. According to the solution, the true anomaly is υ=-0.65° compared to υ=-9.11° of Apollo 12 landing plan. The optimized initial decent position and velocity on the specified orbit is described in Table 8. Duration of this OPLO optimization problem with free descent position on the orbit is tf =477.04s. And the fuel consumption is 6912.91kg. The results indicate that the additional freedom leads to reduction of fuel consumption by ∆m =49.52kg and the flight time by ∆t =160.11s. Table 8 Optimized Initial Decent State on Specific Orbit (OPLO Mission) Variable h0 Vx,0

Value 15125.74m -0.51m/s

Variable α0 Vy,0

Value -9.73° -1645.08m/s

Variable β0 Vz,0

Value -6.48° 406.45m/s

0 -5

Longitude ,deg

Altitude H,m

15000

10000

5000

-10 -15 -20 -25

0 0

200

400 600 Simulation Time,s (a)Altitude H

-30

800

0

200

400 Simulation Time,s (b)Longitude 

600

800

0

200

400 Simulation Time,s (d)Velocity Vx

600

800

0

200

400 Simulation Time,s (f)Velocity Vz

600

800

0

200

400 Simulation Time,s (h)Yaw 

600

800

0 0

Velocity Vx,m/s

Latitude ,deg

-2 -4 -6

-40

-60

-8

-80

-10

0

200

400 Simulation Time,s (c)Latitude 

600

800

0

400

-500

Velocity Vz,m/s

Velocity Vy,m/s

-20

-1000

300 200 100

-1500 0 0

200

400 600 Simulation Time,s (e)Velocity Vy

800

100

0

80

-5

60

-10

Yaw  ,deg

Pitch ,deg

-2000

40 20 0 -20

-15 -20 -25

0

200

400 Simulation Time,s (g)Pitch 

600

800

-30

4

5

x 10

Thrust T,N

4 3 2 1 0 0

200

400 Simulation Time,s (i)Thrust T

600

800

Fig. 7 Optimized Landing Trajectories (OPLO Mission)

7.3. Constrained VTVL with Obstacle Avoidance When the lunar probe performs long-distance transfer maneuvering missions, the flight area is often a complex mountainous terrain. Using the formula (24) to simulate the complex mountainous terrain, the terrain parameters are shown in Table 9. Besides, the total number of Lunar mountains is six (M=6). The margin of safety of ∆hmargin is 200m. Height adjustment parameter of ∆h is 10m. All of the above parameters are given randomly, which are consistent with the detected real lunar terrain. Table 9 Simulation Parameters of Complex Mountainous Terrain

i 1 2 3 4 5 6

hi 5.5km 5.5km 4.5km 4km 5.3km 3km

α0i -23.52° -23.26° -23.4° -23.05° -23.42° -23.7°

β0i -2.72° -2.71° -2.4° -2.3° -2.15° -2.0°

αti 1.52  107 9.1  106 3.64  107 1.52  107 1.52  107 6.06  106

βti 6.06  6.06  1.52  6.06  6.06  6.06 

106 106 107 106 106 106

This scenario uses the final state values of the OPLO as the initial conditions, which are listed in Table 10. The terminal conditions are given in Table 11. β0 and βf in these two tables indicate the lander flies toward the north at a latitude of 1° in the mission. And the attitudes of the lander are vertical at both t0 and tf, which signify the pitch as ϕ0 =90° and ϕf =90°. The obstacle avoidance is formulated as in Eq. (25). Table 10 Initial Conditions (VTVL Mission) Variable h0 Vx,0 m0

Value 0m 0m/s 8321.09kg

Variable α0 Vy,0 ϕ0

Value -23.45° 0m/s 90°

Variable β0 Vz,0

Value -2.94° 0m/s

Table 11 Terminal Conditions (VTVL Mission) Variable hf V⊥,f

Value 0m 0m/s

Variable αf V‖,f

Value -23.45° 0m/s

Variable βf ϕf

Value -1.94° 90°

-23.44

7000 6000

-23.45

Longitude ,deg

Altitude H,m

5000 4000 3000 2000 1000

-23.46

-23.47

0 -1000

0

50

100 150 Simulation Time,s (a)Altitude H

200

-23.48

250

0

50

100 150 Simulation Time,s (b)Longitude 

200

250

100

-2.2

50

Velocity Vx,m/s

Latitude ,deg

-2

-2.4 -2.6

0 -50

-2.8 -100

-3 0

50

100 150 Simulation Time,s (c)Latitude 

200

250

100 150 Simulation Time,s (d)Velocity Vx

200

250

0

50

100 150 Simulation Time,s (f)Velocity Vz

200

250

0

50

150 Velocity Vz,m/s

10

Velocity Vy,m/s

50

200

20

0

-10

-20

0

100 50 0

0

50

100 150 Simulation Time,s (e)Velocity Vy

200

-50

250

100 100

90

50 Yaw  ,deg

Pitch ,deg

80 70 60

-50

50 40 30

0

-100 0

50

100 150 Simulation Time,s (g)Pitch 

200

250

100 150 Simulation Time,s (h)Yaw 

200

250

4

5

x 10

Thrust T,N

4 3 2 1 0 0

50

100 150 Simulation Time,s (i)Thrust T

200

250

Fig. 8 Optimized Landing Trajectories (VTVL mission)

B

B A

A

B

A B A

Fig. 9 Three-Dimensional optimized Trajectory of VTVL Mission

In the VTVL missions, including the unconstrained VTVL mission (VTVL without terrain constraint) and constrained VTVL mission (VTVL with terrain constraints), the iterative calculation procedures are almost same with the steps in the TPLP and OPLP missions, It fails to obtain solution of the two VTVL missions without the initial guess strategy of constraints enhancement. In addition, it takes 28.76 seconds for the unconstrained VTVL mission and 102.37 seconds for the constrained VTVL mission without the initial guess strategy of collocation points refinement. However, when the unified trajectory optimization framework as well as the application of initial guessing strategies are used to calculate the VTVL missions, it takes 12.25 seconds for the unconstrained VTVL mission. With the results of the unconstrained VTVL mission as the initial guess values, the calculation time is 63.52 seconds in total for the constrained VTVL mission. The optimized trajectories of VTVL with no terrain constraint are shown in Fig. 8 with black curves. The duration of VTVL mission is tf =237.14s and the fuel consumption is 1324.63kg. Profiles of the VTVL with terrain obstacle constraints are given in Fig. 8 by blue curves. In order to avoid the obstacle, the lander reaches the height more than 6.5 km (Fig. 8 (a)) before descents. When lander landing, the change in the yaw is approximately ∆ψ=180 ° (Fig. 8 (h)). And the verticality requirements at the initial and final time are reflected by the pitch angle (Fig. 8 (g)). Additionally, both of the two thrust profiles in Fig. 8 (i) satisfiy the discontinuous magnitude constraints Eq. (20). Along the optimal

paths, the thrust switches from the two extremes (minimal or maximal) values in most cases and sometimes the thrust is throttleable between minimal thrust and 65% of full thrust. The whole landing mission takes 257.43 seconds and the fuel consumption is 1349.79kg. It is obvious that with obstacle avoidance the lander takes more flight time and consumes more fuel. The 3D views of the optimal trajectory are given by Fig. 9 in different perspective. In this figure, points A and B are the initial position and final destination, respectively. The constrained VTVL trajectories with obstacle avoidance are represented by blue curves, and the nominal VTVL trajectories are represented by black curves. It is clear that the blue trajectories avoid several mountains and securely land at the destination. 8. Conclusion In this study, minimum-fuel trajectories have been optimized for various landing missions with specific requirements and complex constraints. A unified trajectory optimization framework based on GPM is proposed for solving the underlying optimal control problems numerically. Advanced initialization strategies are developed to improve the performance of optimization when dealing with complex constraints, such as discontinuous thrust, obstacle avoidance, multi-phase connection, and boundary constraints. Three scenarios of lunar soft landing have been investigated and the trajectories are compared. For the TPLP mission, OPLP mission, as well as the VTVL missions, the results demonstrate the performance and the adaptability of the proposed trajectory optimization framework in solving difficult problems derived from various landing missions. The proposed unified framework is combined with the application of advanced initialization strategies to generate the initial guess values of optimization proposition, which can guarantee the success of various lunar soft landing missions, enhance the robustness of convergence, and improve the efficiency of the optimization calculation. Moving-horizon closed-loop optimization accounting for real-time disturbances and uncertainties would be investigated for the future work. In addition, the unified trajectory optimization framework for lunar soft landing is open for other complex features and constraints that are not designed in this paper. It is flexible enough to deal with different complex constraints that are described as equality or inequality constraints. Recently, spacecraft systems are required to have autonomy and support autonomous decision-making and intelligent control, which are very important for the guidance and control of spacecraft. The unified trajectory optimization framework for lunar soft landing is beneficial to the research and understanding in the early stage, such as the proposition construction, model establishment and problem analysis. And eventually it can promote successful implementation of the lunar surface autonomous soft landing missions. Acknowledgments This research was supported by the 973 Program of China (No. 2012CB720503) and the National Nature Science Foundation of China (No. 61203132).

References NASA Staff, “Apollo Program Summary Report,” NASA Rept. JSC-09423, Apr. 1975. Boere M. Optimization of descent trajectories for lunar base settlement. Faculty of Aerospace Engineering Astrodynamics & Space Missions, Delft University of Technology; 2010. Master Dissertation. [3] Young A. The twenty-first century commercial space imperative. New York: Springer; 2015. p. 43–57. doi:10.1007/978-3-319-18929-1. [4] F. V. Bennett, T. G. Price, “Study of Powered-Descent Trajectories for Manned Lunar Landings,” NASA TN D-2426, Aug. 1964. [5] D. C. Cheatham, F. V. Bennett, T. M. Branch, “Apollo Lunar Module Landing Strategy,” Apollo Lunar Landing Mission Symposium, NASA TM X-58006, pp. 175-240, 1966. [6] G. Elverum, P. Staudhammer, J. Miller, A. Hoffman, R. Rockow, “The Descent Engine for the Lunar Module”, 3rd Propulsion Joint Specialist Conference, No. 67-521, 1967. [7] G. A. Dressler, “Summary of Deep Throttling Rocket Engines with Emphasis on Apollo LMDE,” 42nd IAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA 2006-5220, 2006. [8] M. J. Policelli, “Vertical Takeoff Vertical Landing Spacecraft Trajectory Optimization via Direct Collocation and Nonlinear Programming”, Master of Science Thesis, Department of Aerospace Engineering, the Pennsylvania State University, Aug. 2014. [9] Nigro N. Knack night sky: decoding the solar system, from constellations to black holes. Globe Pequot Press; 2010. p. 80. [10] J. T. Betts, “Survey of Numerical Methods for Trajectory Optimization,” Journal of Guidance, Control, and Dynamics, Vol. 21, No. 2, pp. 193-207, Mar. 1998. [11] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelize, and, E. F. Mishchenko, “The Mathematical Theory of Optimal Processes,” Wiley, New York, 1962, Chap. 2. [12] D. E. Kirk, “Optimal Control Theory: An Introduction,” Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1970. [1] [2]

[13] J. T. Betts, “Practical Methods for Optimal Control Using Nonlinear Programming” 2nd ed., Society for Industrial and Applied Mathematics, Philadelphia, 2010. [14] D. A. Benson, “A Gauss Pseudospectral Transcription for Optimal Control,” Ph.D. Dissertation, Dept. of Aeronautics and Astronautics, Massachusetts Inst. of Technology, Cambridge, Feb. 2005. [15] D. A. Benson, G. T. Huntington, T. P. Thorvaldsen, and A. V. Rao, “Direct Trajectory Optimization and Costate Estimation via an Orthogonal Collocation Method,” Journal of Guidance, Control, and Dynamics, Vol. 29, No. 6, pp. 1435–1440, Nov. 2006. [16] G. T. Huntington, “Advancement and Analysis of a Gauss Pseudospectral Transcription for Optimal Control Problems,” Ph.D. Dissertation, Dept. of Aeronautics and Astronautics, Massachusetts Inst. of Technology, Cambridge, MA, May 2007. [17] Fahroo F, Ross IM. Pseudospectral methods for infinite-horizon nonlinear optimal control problems. J Guid Control Dyn 2008;31(4):927–36. doi:10.2514/1.33117. [18] Garg D, et al. Direct trajectory optimization and costate estimation of finitehorizon and infinite-horizon optimal control problems using a radau pseudospectral method. Comput Optim Appl 2009;49(2):335–58. doi:10.1007/s10589-009-9291-0. [19] Kameswaran S, Biegler LT. Convergence rates for direct transcription of optimal control problems using collocation at radau points. Comput Optim Appl 2008;41(1):81–126. doi:10.1007/s10589-007-9098-9. [20] Garg D, et al. A unified framework for the numerical solution of optimal control problems using pseudospectral methods. Automatica 2010;46(11):1843–51.doi:10.1016/j.automatica.2010.06.048. [21] G. T. Huntington, A. V. Rao, “Optimal Reconfiguration of Spacecraft Formations Using the Gauss Pseudospectral Method,” Journal of Guidance, Control, and Dynamics, Vol. 31, No. 6, pp. 689–698, 2008. [22] T. R. Jorris, C. S. Schulz, F. R. Friedl, A. V. Rao, “Constrained Trajectory Optimization Using Pseudospectral Methods,” AIAA Atmospheric Flight Mechanics Conference and Exhibit, 2008. [23] T. R. Jorris, R. G. Cobb, “Three-dimensional Trajectory Optimization Satisfying Waypoint and No-fly Zone Constraints,” Journal of Guidance, Control, and Dynamics, Vol. 32, No. 2, pp. 551-572, 2009. [24] B. Tian, W. Fan, R. Su, Q. Zong, “Real-time trajectory and attitude coordination control for reusable launch vehicle in reentry phase,” IEEE Trans. Ind. Electron., vol. 62, no. 3, pp. 1639-1650, 2015. [25] F. Zhang, G. Duan, “Integrated translational and rotational control for the terminal landing phase of a lunar module,” Aerospace Science and Technology, Vol. 27, No. 1, pp. 112-126, 2013. [26] M. Pontani, G. Cecchetti, P. Teofilatto, “Variable-Time-Domain Neighboring Optimal Guidance, Part 2: Application to Lunar Descent and Soft Landing,” Journal of Optimization Theory and Applications, Vol. 166, No. 1, pp. 93-114, 2015. [27] B. Zhang, S. Tang, B. Pan, “Multi-constrained suboptimal powered descent guidance for lunar pinpoint soft landing,” Aerospace Science and Technology, Vol. 48, pp. 203-213, 2016. [28] V. L. Rvachev, T. I. Sheiko, “R-Functions in Boundary Value Problems in Mechanics,” Applied Mechanics Reviews, Vol.48, No.4, pp.151-188, 1995. [29] L. Kurpa, V. Rvachev, E. Ventsel, “The R-Function Method for the Free Vibration Analysis of Thin Orthotropic Plates of Arbitrary Shape,” Journal of Sound and Vibration, Vol. 261, No. 1, pp. 109-122, 2003. [30] N. X. Vinh, A. Busemann, R. D. Culp, “Hypersonic and Planetary Entry Flight Mechanics,” Ann Arbor, MI:University of Michigan Press, 1980． [31] O. Montenbruck, E. Gill, “Satellite orbits: models, methods and applications,” Springer Science & Business Media, 2012. [32] P. E. Gill, W. Murray, and M. A. Saunders, “SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization,” SIAM Journal on Optimization, Vol. 12, No. 4, pp. 979–1006, 2002. [33] G. L. Carman, M. N. Montez, “Apollo 12 (Mission H-1) Spacecraft Dispersion Analysis Volume Ⅳ: Descent and Ascent Dispersion Analyses Part 1: Lunar Descent,” MSC Internal Note No. 69-FM-275, NASA-TM-X-72128, Manned Spacecraft Center, Oct. 1969. [34] W. Joe, Nolley, “Apollo 12 (Mission H-1) Spacecraft Dispersion Analysis Volume Ⅶ: Dispersion Summary,” MSC Internal Note No. 69-FM-292, NASA-TM-X-72099, Manned Spacecraft Center, Nov. 1969. [35] F. V. Bennett, “Apollo Experience Report-Mission Planning for Lunar Module Descent and Ascent,” NASA Technical Note, NASA TN D-6846, Manned Spacecraft Center, Jun. 1972.

   

A unified trajectory optimization framework is proposed for lunar soft landing. Advanced initialization strategies improve the performance of optimization. Use R functions to describe discontinuous thrust constraints in the same problem. Results demonstrate the adaptability of the framework for various landing missions.