- Email: [email protected]

Corso di Laurea in Ingegneria Spaziale

ROBUST CONTROL FOR PLANETARY LANDING MANEUVERS

Relatore: Prof. Michèle LAVAGNA Correlatore: Ing. Roberto ARMELLIN

Tesi di laurea di: Paolo LUNGHI Matr. 734407

Anno Accademico 2011 - 2012

A Papà.

Abstract This work focuses on an adaptive guidance algorithm for planetary landing that updates the trajectory to the surface by means of a minimum fuel optimal control problem solving. A semi-analytical approach is proposed. The trajectory is expressed in a polynomial form of minimum order to satisfy a set of 17 boundary constraints: 12 constraints on initial and ﬁnal state and 5 control constraints, added in order to include attitude requirements. By imposing boundary conditions, a fully determined guidance proﬁle is obtained, function of only two parameters: time-of-ﬂight and initial thrust magnitude. The optimal guidance computation is reduced to the determination of these parameters, according to additional path constraints due to the actual lander architecture: available thrust and control torques, visibility of the landing site, and other additional constraint not implicitly satisﬁed by the polynomial formulation. Solution is achieved with a simple two-stage compass search algorithm: the algorithm ﬁrstly ﬁnds a feasible solution; whenever detected, it keeps solving for the optimum; nonlinear constraints are evaluated numerically, by pseudospectral methods. Results on diﬀerent scenarios for a Moon landing mission are shown and discussed to highlight the eﬀectiveness of the proposed algorithm and its sensitivity to the navigation errors. Keywords: pinpoint landing, adaptive guidance, retargeting, hazard avoidance

Sommario Il presente lavoro è focalizzato sulla formulazione e la veriﬁca di un algoritmo di guida adattiva per l’atterraggio planetario, che in seguito alla modiﬁca del sito di atterraggio riformuli la traiettoria risolvendo un problema di ottimizzazione del carburante. Viene qui proposto un approccio semi-analitico. La traiettoria viene rappresentata con una forma polinomiale, del minimo ordine necessario per soddisfare un set di 17 condizioni al contorno, di cui 12 per gli stati iniziale e ﬁnale del sistema e 5 riguardanti le variabili di controllo, imposte dall’assetto iniziale e ﬁnale desiderato. Imponendo queste condizioni al contorno, separatamente su ogni asse, si ottiene un completo proﬁlo di guida, funzione di soli due parametri, tempo di volo e spinta iniziale. Il problema è quindi ridotto all’individuazione dei valori di questi parametri, tali da ottimizzare il consumo di propellente soddisfacendo al tempo stesso i vincoli addizionali posti dall’architettura del lander: valori minimi e massimi erogabili di spinta e coppie di controllo, requisiti di visibilità sul sito di atterraggio e tutti gli eventuali altri vincoli non implicitamente soddisfatti dalla formulazione polinomiale. La soluzione del problema di ottimo viene ottenuta mediante un semplice algoritmo di compass search in due stadi: da principio l’algoritmo procede alla ricerca di una soluzione che non violi i vincoli; una volta individuata, prosegue alla ricerca dell’ottimo. I vincoli non lineari vengono valutati discretamente, mediante metodi pseudospettrali. L’algoritmo viene quindi applicato a diﬀerenti scenari di un atterraggio lunare, in modo da stimarne ﬂessibilità di applicazione e sensitività agli errori di navigazione. Parole chiave: atterraggio di precisione, guida adattiva, retargeting, hazard avoidance

Ringraziamenti

Un ringraziamento alla professoressa Lavagna, che mi ha dato l’opportunità di svolgere questa tesi e nei cui corsi sono rimasto coinvolto nel problema dell’atterraggio. Grazie all’Ing. Roberto Armellin che, con i suoi consigli al momento giusto, è stato determinante per la buona riuscita del lavoro. Grazie alla mia famiglia, a mia madre, alle mie sorelle Mari, Claudia e Daniela, a Simone, Angelo e Alberto, che mi hanno sempre sostenuto in ogni momento e in ogni situazione, e sono sempre stati per me il migliore degli esempi. Grazie a tutti i miei nipotini, Michele, Milena, Andrea, Giulia, Arianna, Marta e Edoardo: distraendo lo zio Paolo che lavorava, gli hanno consentito di riposare i neuroni e ripartire ogni volta più scattante di prima. Grazie al Ricky, al Dimpo, al Manu, a Stephane a Michel e al Mike, senza i quali l’università sarebbe stata davvero più pesante. Grazie al Rosso, al Toga e al Betto, che han sempre sopportato le mie sparizioni prolungate. Grazie infine a Pao, pungolo costante a dare il meglio di me e che in me ha sempre creduto.

Contents

Nomenclature

xvii

1 Introduction 1.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Motivations and Goal . . . . . . . . . . . . . . . . . . . . . . . 1.3 Dissertation overview . . . . . . . . . . . . . . . . . . . . . . . 2 Nominal Landing Trajectory 2.1 Landing phases . . . . . . . . . . . . . . . . 2.2 Nominal Landing Simulation . . . . . . . . . 2.2.1 Discretization of Control Variables . 2.2.2 Nominal Landing Simulation Results 3 Retargeting Algorithm 3.1 Landing Model . . . . . . . . . 3.1.1 Reference systems . . . . 3.1.2 Equation of motions . . 3.2 Trajectory Design . . . . . . . . 3.2.1 Polynomial Formulation 3.2.2 Additional Constraints . 3.3 Optimization Algorithm . . . . 3.4 Algorithm Performances . . . . 3.4.1 Order of approximation 3.4.2 Divert Capability . . . . 3.4.3 Optimality . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

1 2 3 4

. . . .

7 7 9 11 15

. . . . . . . . . . .

19 19 19 20 23 24 27 30 34 35 37 39

4 Retargeting Simulation 45 4.1 Modeling criteria . . . . . . . . . . . . . . . . . . . . . . . . . 46 xi

4.2

4.1.1 Lander Dynamics . . . . . . . . 4.1.2 Navigation Errors Model . . . . 4.1.3 Guidance . . . . . . . . . . . . 4.1.4 Control and Actuation . . . . . Results . . . . . . . . . . . . . . . . . . 4.2.1 Single Landing Case . . . . . . 4.2.2 Sensitivity to TLS . . . . . . . 4.2.3 Sensitivity to Initial Dispersion 4.2.4 Cameras Field of View . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

46 50 52 54 58 58 63 65 67

5 Conclusions 71 5.1 Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2 Future Developments . . . . . . . . . . . . . . . . . . . . . . . 72 Bibliography

xii

75

List of Figures

1.1 Historical Perspective over the Landing Accuracy on Mars . .

2

2.1 2.2 2.3 2.4 2.5 2.6 2.7

Landing phases. . . . . . . . . . . . . . . . . . . . . . Global landing reference system. . . . . . . . . . . . . Nominal landing - control proﬁles . . . . . . . . . . . Nominal landing path, Vertical and horizontal speed . Nominal landing path, altitude and mass . . . . . . . Nominal landing - view angle onto NLS . . . . . . . . Nominal landing - Trajectory . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

8 9 15 16 17 17 18

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16

Body-Fixed, Flight and Ground reference systems. Body-Fixed reference system and euler angles. . . Glide-slope constraint. . . . . . . . . . . . . . . . Considered RCS thrusters scheme. . . . . . . . . . ESA Lunar Lander dimensions. . . . . . . . . . . Mean Computation time vs N . . . . . . . . . . . Position precision vs N . . . . . . . . . . . . . . . Speed precision vs N . . . . . . . . . . . . . . . . Feasibility performances. . . . . . . . . . . . . . . Feasibility performances (2) . . . . . . . . . . . . Computational performances. . . . . . . . . . . . Fuel consumption performances. . . . . . . . . . . Optimality comparison, CASE 1. . . . . . . . . . Optimality comparison, CASE 2. . . . . . . . . . Optimality comparison, CASE 3. . . . . . . . . . Optimality comparison, CASE 3. . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

20 21 30 34 35 36 37 37 38 39 39 40 41 42 43 44

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

4.1 Complete landing system logical scheme . . . . . . . . . . . . 45 4.2 Simpliﬁed Landing System . . . . . . . . . . . . . . . . . . . . 46 xiii

4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 4.19

xiv

Parametric lander inertial model . . . . . . . . . . . . . Guidance logic diagram . . . . . . . . . . . . . . . . . . Control logic diagram . . . . . . . . . . . . . . . . . . . PWPF Modulation . . . . . . . . . . . . . . . . . . . . Single Retargeting, Position . . . . . . . . . . . . . . . Single Retargeting, Speed . . . . . . . . . . . . . . . . Single Retargeting, target and actual attitude . . . . . Single retargeting. Angular velocities . . . . . . . . . . Single retargeting. Mass trend . . . . . . . . . . . . . . Single retargeting. ACS Thrusters ﬁrings . . . . . . . . Single retargeting. Main thrust magnitude trend . . . . Landing Position Sensitivity to the TLS . . . . . . . . Landing Velocity Sensitivity to the TLS . . . . . . . . Landing Position Sensitivity to Navigation dispersions . Landing Velocity Sensitivity to Navigation dispersions . Camera pointing during retargeting . . . . . . . . . . . Sightline-TLS angle of view . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

48 53 54 57 58 59 60 61 61 62 62 63 64 65 66 67 68

List of Tables

2.1 Nominal landing simulation assumptions. . . . . . . . . . . . . 15 2.2 Nominal landing simulation - Additional constraints values. . . 16 3.1 Lander architecture assumptions. . . . . . . . . . . . . . . . . 35 4.1 4.2 4.3 4.4 4.5 4.6

Components of lander inertial model. . . . . . Sensors performance properties. . . . . . . . . Navigation errors covariance at 2000 m. . . . . Navigation errors covariance at ground (0 m). PID gains and equivalent dynamics properties. PWPF Modulator parameters. . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

48 51 52 52 56 57

xv

Nomenclature

∆

Compass search mesh size

Φ

Compass search Modiﬁed cost function

β

Angle of thrust

β

Angle of thrust pseudospectral representation vector

φ, θ, ψ Euler angles: Roll, Pitch and Yaw γ

Glide-slope angle

µ

Moon standard gravitational parameter

ω

Angular speed vector

ωM

Measured angular speed vector

ωT

Target angular speed vector

ωe

Error angular speed vector

ψ

Yaw angle pseudospectral representation vector

ρ

Margin of safety weight on maximum control torques

τ

Pseudospectral computational time

τ∆

Compass search minimum mesh size

τk

Chebyshev-Gauss-Lobatto nodes

τm

PWPFM ﬁlter time constant xvii

θM

Azimuth in polar coordinates, Moon centered frame

θ

Pitch angle pseudospectral representation vector

ADC

Attitude director cosine matrix

AGF

Director cosine matrix coordinates transformations between Flight and Ground frames

B

Gyro bias error

DN

Chebyshev diﬀerentiation matrix

˜N D

Modiﬁed Chebyshev diﬀerentiation matrix

F

Feasibility function

gM

Gravitational acceleration on Moon’s surface

I

Matrix of inertia

Imax

Maximum moment of inertia

Isp

Speciﬁc impulse

Kd

PID derivative gain

Ki

PID integral gain

Km

PWPFM ﬁlter gain

Kp

PID proportional gain

M

Gyro misalignment error

MC

Control torque

MC

Control torques vector

˜C M

Theoretical control torques vector

MD

Disturbance torques vector

Mmax Maximum available control torque N

Order of pseudospectral approximation

P

Thrust-to-mass ratio

xviii

P

Cartesian Thrust-to-mass ratio vector

RM

Moon mean radius

S

Gyro scale factor error

T

Thrust magnitude

T

Cartesian thrust vector

TT

Target thrust magnitude

Tmag Thrust magnitude pseudospectral representation vector a

Semi-major axis of orbit

a

Cartesian acceleration vector

b

Thrust point of application vector

cl , cU Non-linear constraints Lower and Upper bounds c

Optimization non-linear constraints vector

˜ c

Generalized non-linear constraints vector

ea

Euler angles vector

f

Optimization cost function

g0

Standard gravity acceleration

g

Cartesian gravity acceleration vector

h

Altitude

itF max Maximum feasibility iterations limit itOmax Maximum optimality iterations limit m

Mass of the spacecraft

mM

Measured spacecraft mass

mPDI Spacecraft mass at Powered Descent Initiation mdry Dry mass of the spacecraft mfuel Total fuel mass on the spacecraft xix

n

Gyro sensor noise

p, p

Generic function of time, and vector of its discrete pseudospectral representation

q

Quaternions vector

qM

Measured quaternions vector

qT

Target quaternions vector

qe

Error quaternions vector

r

Radial distance from origin in polar coordinates, Moon centered frame

r

Cartesian position vector

rM

Measured position vector

t

Time

t0 , tf Initial and ﬁnal time tToF

Time-of-ﬂight

u

Tangential speed in polar coordinates, Moon centred frame

umax PWPFM Schmitt-trigger output uoff

PWPFM Schmitt-trigger oﬀ-value

uon

PWPFM Schmitt-trigger on-value

v

Radial speed in polar coordinates, Moon centered frame

v

Cartesian speed vector

vM

Measured speed vector

wC

Clenshaw-Curtis quadrature scheme weights

wF

Generalized non-linear constraints weights

xL , xU Lower and Upper bounds of optimization variables x

Optimization variables vector

˜ x

Normalized optimization variables vector

xx

B

Body-ﬁxed reference frame

F

Flight axes reference frame

G

Ground axes reference frame

ACS Attitude Control System AG

Approach Gate

CGL Cebyshev-Gauss-Lobatto CoM Center of Mass DOI Descent Orbit Initiation FoV

Field of View

HDA Hazard Detection and Avoidance HG

High Gate

IMU Inertial Measurement Unit LG

Low Gate

LLO Low Lunar Orbit MB

Main Brake

NLS Nominal Landing Site PDI

Powered Descent Initiation

TD

Terminal Descent

TLS Target Landing Site

xxi

Chapter

1

Introduction In last years, a renewed interest in planetary exploration has brought to the realization of several missions, especially towards Mars, culminated with the recent landing of the rover Curiosity in August 2012. Together with Mars, the Moon is a main destination for exploration. The European Space Agency has conducted several studies concerning a possible unmanned lunar lander, while NASA is planning to send humans back to the Moon. ESA will supply the Orion/MPCV Service Module (SM) for the 2017 unmanned Exploration1 Mission, including ground and ﬂight operation support. Provisions for the construction and delivery of a second SM have been taken. Recently ESA and the Russian federal space agency, Roscosmos, have signed a formal agreement to work in partnership on the ExoMars programme towards the launch of two missions in 2016 and 2018, in order to bring a rover on Mars surface. In all these missions, the Entry Descent and Landing phase fulﬁls a critical role. During last decades, several improvements in automatic landing precision have been done (see Figure 1.1), but the relative uncertainty over the ﬁnal achieved position still imposes strict requirements onto the landing site choice. This is why an autonomous, precise and safe landing capability is a key feature for the next space systems generation. A dynamical landing site selection could allow the reaching of more scientiﬁcally relevant places, avoiding eventual hazardous terrain items. The short duration of the terminal descent phase and the delay of communications at large distances make a direct ground control impossible, and impose the develop of a fully autonomous system. After the landing site selection, the system needs to recalculate a pinpoint feasible trajectory toward the target. The pinpoint landing problem can be deﬁned as guiding a lander spacecraft to a given target on the planet’s surface 1

Chapter 1 with an accuracy of fewer than several hundred meters [2]. This implies the resolution of a minimum fuel optimal control problem.

Figure 1.1: This image illustrates how spacecraft landings on Mars have become more and more precise over the years. Since NASA’s first Mars landing of Viking in 1976, the targeted landing regions, or ellipses, have shrunk. Image credit: NASA/JPL-Caltech/ESA.

1.1

State of the art

A trajectory based on a quartic polynomial in time was used during the Apollo missions [19]. This allowed to perform in 1969 the ﬁrst recognized pinpoint landing (also if manually aided) during the Apollo 12 mission [4], whose lunar module landed at only 600 m from the target, the Surveyor 3 probe landing site. A derivative of the Apollo lunar descent guidance was still considered in recent years, for the Mars Science Laboratory (MSL) [35]. Various other approaches to obtain both numerical and approximate solutions of the pinpoint landing terminal guidance problem have been described over the last few years. In [32] the ﬁrst-order necessary conditions for the problem are developed, and it is shown that the optimal thrust proﬁle has a maximum-minimum-maximum structure. 2

Introduction Direct numerical methods for trajectory optimization have been widely investigated, not requiring the explicit consideration of the necessary conditions [5]. These methods have been used together with Chebyshev pseudospectral techniques, in order to reduce the number of the optimization variables [10]. Also convex programming approach has been proposed, in order to guarantee the convergence of the optimization [2]. Direct collocation methods has showed that the size of the region of feasible initial states, for which there exist feasible trajectories, can be increased drastically (more than twice) compared to the traditional polynomial-based guidance approaches, but at the price of a higher computational cost [2].

1.2

Motivations and Goal

The aim of this work is the development of a guidance algorithm capable to dynamically recompute and correct the landing trajectory during the descent, allowing the on-board choice of the landing site, required by systems that have to operate in full autonomy. Some principles has been considered as guidelines: • Computational Efficiency A full autonomous retargeting requires algorithms executable on-board, with limited hardware capabilities, and minimal time lag (if not in realtime). Low computational cost is a primary requirement. • Flexibility The attainable landing area must be maximized, in order to have the highest probability to ﬁnd a safe landing site. • Landing Accuracy Landing retargeting and hazard avoidance imply the requirement of a precision pinpoint landing. A semi-analytical approach is proposed, in order to conjugate the low computational cost of polynomial approximation to the larger ﬂexibility of direct optimization methods. Fuel consumption has been used as optimality criterion. In this work the last phases of a lunar landing have been examined. This implies that no aerodynamic forces are considered. However, the obtained results remain still valid also in presence of atmosphere with low density, such as the case of Mars, due to the relatively small velocities involved in this phase. 3

Chapter 1 In all the tests performed the ESA Lunar Lander mission has been taken as reference scenario. The European Lunar Lander is a mission under study within the Human Spaceﬂight and Operations Directorate of the European Space Agency (ESA). Originally planned for launch in 2018 and designed for landing near the Moon’s south pole, the mission’s primary objectives include the demonstration of safe precision landing technology as part of preparations for participation to future human exploration of the Moon [26]. Recently, the project was put on hold at the 2012 ESA Ministerial Council [8]. The technology developed in the context of Lunar Lander phase B1 could be exploited for future cooperations in the area of Lunar Exploration with Russia. The Luna-Resource Lander mission, planned by Roscosmos for 2017, could be a testing platform for European precision landing technology, with the proposed Hazard Detection and Avoidance Experiment and the Visual Absolute/Relative Terrain Navigation Experiment (VNE). The obtained results could be employed to perform an autonomous precision landing in the Lunar Polar Sample Return (LPSR) mission, planned for launch in 2020 [12].

1.3

Dissertation overview

In the ﬁrst part of Chapter 2 diﬀerent phases of a lunar landing, from parking orbit to touchdown, are presented. In the second part, a fuel-optimal solution for a nominal lunar landing trajectory (needed as starting point for retargeting) is obtained through numerical optimization in a two-dimensional frame. Chapter 3 focuses on the retargeting guidance algorithm development. The ﬁrst part is dedicated to illustrate the reference systems considered. The ODE system that describes the lander translational dynamics is derived. Then, a polynomial approach is used to reduce the trajectory computation to the search of optimal values for only two variables: the time of ﬂight and the initial thrust magnitude. A compass search algorithm, modiﬁed to handle nonlinear constraints, is employed in order to ﬁnd the fuel optimal solution of the problem, with a light computational cost. The last part of the chapter is dedicated to the validation of the algorithm: computational performances and divert capabilities at diﬀerent altitudes are estimated by Monte Carlo analysis. Finally, the optimality of the achieved solution is appraised through comparisons with direct collocation schemes. In Chapter 4 a closed-loop simulation of a retargeting during a lunar landing is presented. First, the spacecraft dynamic system, completed by 4

Introduction the rotational dynamics, is presented. In the second part of the chapter, the model structure is expounded; assumptions and approximations made are discussed. The guidance algorithm developed in Chapter 3 is coupled with an all-thruster Attitude Control System. A vision-based navigation system is emulated by introducing stochastic errors in the states inputs of the guidance and control modules. In the last part, obtained results are illustrated. First, a typical case is presented, in order to show distinctive behaviours in the system response. Then, the landing accuracy is evaluated using Monte Carlo analysis. Sensitivity to the diversion magnitude and to navigation dispersion are investigated. Possible sources of dispersion, and their relative impact over the landing accuracy are discussed. Finally, preliminary considerations about the obtained visibility of the landing site are presented. In Chapter 5 obtained results are summarized, and some considerations on requirements imposed by the retargeting strategy are given. An outline of possible future developments and further investigations closes the work.

5

Chapter

2

Nominal Landing Trajectory In this chapter all the phases of a lunar landing from a LLO parking orbit is presented. Then, a nominal two-dimensional landing trajectory, needed as starting point for the development of a retargeting algorithm, is obtained through a fuel optimization. Because fuel mass is a major part of the total vehicle mass, minimizing the fuel consumption is the most suitable criterion in determining an eﬃcient landing strategy.

2.1

Landing phases

A nominal Descent and Landing proﬁle, as described in ref. [15] is considered. All the landing phases are depicted in Figure 2.1. The lander is assumed to be delivered on a 100 km Low Lunar Orbit (LLO). The landing procedure is made by two main phases, the Coasting and the Powered Descent. The Powered Descent is, in turn, divided in three subphases: the Main Brake, the Approach and the Terminal Descent. • Coasting Phase At Descent Orbit Initiation (DOI), a de-orbit burn inserts the lander in a Descent Orbit, with a periselenium at 15 km altitude. In this phase, the spacecraft needs to execute a 180◦ turn around the pitch axis, in order to point the main engine nozzle at forward direction at the periselenium. • Main Brake Phase At periselenium, the Powered Descent phase is initiated (PDI), starting with the Main Braking phase: the engines are ignited to provide maximum thrust and to eﬃciently remove as much of the orbit velocity 7

Chapter 2 as possible. At the High Gate (HG) point, the Nominal Landing Site (NLS) comes into the ﬁeld of view of the sensors. The lander starts a pitch maneuver to realize a trajectory compatible with navigation and Hazard Detection and Avoidance (HDA) constraints, in terms of visibility of the NLS and time for image processing. • Approach Phase A few seconds later, at Approach Gate (AG) the Hazard Detection System comes into operation. Once the hazard map of the NLS area is completed, the thrust level is reduced to gain manoeuvrability and between 2000 and 100 m of altitude one or two retargetings may be commanded to divert to a safer Landing Site (Target Landing Site, TLS). • Terminal Descent At Low Gate (LG), the Terminal Descent phase begins. This part starts about 30 m above the surface and achieves a vertical uniform trajectory until Touchdown (TD). The thrust is basically equal to the weight of the lander during this last phase.

Figure 2.1: Lunar landing phases with representative timeframe and altitude values (adopted from ref. [15]).

8

Nominal Landing Trajectory

2.2

Nominal Landing Simulation

The landing trajectory is obtained by ﬁnding a fuel-optimal solution of the Main Brake phase and of the subsequent Approach phase, from PDI to Low Gate, towards the Nominal Landing Site. Searching for a nominal trajectory, no retargeting is for now considered. Every position on a planet can be reached with a planar trajectory, if the parking orbit inclination is equal to the ﬁnal desired latitude. Then, the landing is modelled as two-dimensional. Assuming r and θM as polar coordinates of the position of the lander (see Figure 2.2), u as tangential speed, v as radial speed, m as lander mass, the lander translational dynamics is governed by the system: r˙ = v, u θ˙M = , r u˙ = T (t) sin β(t) − 2 uv , m r (2.1) 2 µ u T (t) cos β(t) − 2 + , v˙ = m r r T (t) m ˙ =− , Isp g0

where the thrust magnitude T (t) and the thrust angle β(t) are control variables; µ is the standard gravitational parameter of the Moon, Isp is the speciﬁc impulse of the lander main engine, g0 is the standard gravity acceleration on Earth. Assuming the altitude h(t) = r(t) − RM , where RM is the

u

T β

M oo n s ur f a ce

v

m r ϑ

Figure 2.2: Global landing reference system.

9

Chapter 2 Moon mean radius, indicating with the subscript p the values at the transfer orbit periselenium, the initial system states are: h(0) = hp , θM (0) = 0, (2.2) u(0) = up , v(0) = 0, m(0) = m , p

together with the additional control constraint: π β(0) = − . 2

(2.3)

An initial mass mp of 1500 kg is expected [15]; hp = 15 000 m is the transfer orbit altitude at periselenium, and up is the corresponding orbital velocity: s 2 1 = 1696.02 m s−2 . up = µ − rp a Since a polar landing is contemplated, no horizontal speed is required at the touchdown. The ﬁnal states lower and upper bounds considered are: m, 20 ≤ h(tf ) ≤ 30 −0.5 ≤ u(tf ) ≤ 0.5 m s−1 , −1.5 ≤ v(tf ) ≤ 0 m s−1 ,

(2.4)

along with the ﬁnal control constraint of vertical attitude: β(tf ) = 0.

(2.5)

Final mass m(tf ) and polar angle θM (tf ) are left unconstrained. In order to obtain a fuel-optimal solution, the depleted propellant ffuel = mPDI − m(tf ) is the optimization cost function. In addition with these boundary constraints, the landing proﬁle must satisfy some additional constraints during the descent: • The lander must not hit the ground: h(t) ≥ 0. 10

(2.6)

Nominal Landing Trajectory • Requested control torques must be compatible with the actual lander architecture: −Mmax ≤ MC (t) ≤ Mmax . In a planar problem, the spacecraft attitude is controlled only around the pitch axis. Assuming the absence of thrust vectoring control, the angle of thrust β can be exploited to obtain requested control torque. If I is the moment of inertia around pitch axis: ˙ MC = I β¨ + I˙β. In this simulation, the moment of inertia is assumed constant and equal to Imax , that is its maximum value, at the beginning of the simulation. This leads to an overestimation of requested control torques. The resulting control torque constraint is: − Mmax ≤ Imax β¨ ≤ Mmax .

(2.7)

• The thrust vector must always point away from the ground: − π/2 ≤ β(t) ≤ π/2.

(2.8)

• The requested thrust magnitude have to match the on-board engine capabilities: Tmin ≤ T (t) ≤ Tmax . (2.9) The resulted problem is: Find the suitable control profiles β(t) and T (t) and the maneuver time tf that, with initial conditions (2.2), applied to the dynamic system (2.1), minimize the cost function ffuel , according to constraints (2.4), (2.3), (2.5), (2.6), (2.7), (2.8) and (2.9).

2.2.1

Discretization of Control Variables

In this form, the control proﬁles are continuous variables, with consequently inﬁnite values assumed in each time interval. To make them easier to handle by optimization programs, they have been discretized by pseudospectral techniques. If p(t) is a generic function of time deﬁned in the interval [tini , tend ], it can be represented by a ﬁnite set of values pk , k = 1, 2, . . . , N, that are the values assumed by the function at Chebyshev-Gauss-Lobatto (CGL) nodes 11

Chapter 2 points (see [6, 25]). CGL points lie in the computational interval [−1, 1], and are deﬁned as: τk = − cos(πk/N), k = 0 . . . N. (2.10)

The function evaluation is shifted onto this interval by the linear transformation: t ∈ [tini , tend ] ⇒ τ ∈ [−1, 1], t(τ ) = [(tini − tend )τ + (tini + tend )]/2,

(2.11)

τ ∈ [−1, 1] ⇒ t ∈ [tini , tend ], τ (t) = [2t − (tini + tend )]/(tend − tini ).

(2.12)

Then, the vector p is made by the function values at corresponding points: p0 p1 pk = p t(τk ) . (2.13) p = .. , . pN

The function can be evaluated at each time instant with a polynomial approximation of the form: p(t) ≃ p˜(t) =

N X

pk ϕk (τ (t)),

(2.14)

k=0

where ϕk (τ (t)), for k = 0, 1, . . . , N are the Lagrange interpolating polynomials of order N. ¨ The second derivative β(t), necessary for the evaluation of the control torque constraint, is obtained by the use of the Chebyshev diﬀerentiation matrix. The derivative β˙ in the computational domain at CGL points is obtained through a simple matrix product: β˙ = DN β. The matrix components are: ci (−1)i+j , cj τi − τj −τj , 2(1 − τj2 ) (DN )ij = 2N 2 + 1 , 6 2 − 2N + 1 , 6

12

(2.15)

i 6= j, 1 ≤ i = j ≤ N − 1, i = j = 0, i = j = N.

(2.16)

Nominal Landing Trajectory τj = cos(πj/N) are CGL points of Nth order, and ( 2, cj = 1,

j = 0, N 1 ≤ j ≤ N − 1.

Cebyshev matrix diﬀerentiates a function in the computational domain. To obtain derivatives in the real domain, a change of variable is needed. From eqn. (2.11): tend − tini dτ, 2 ⇓ dβ 2 dβ 2 DN β. = = dt tend − tini dτ tend − tini dt =

(2.17)

(2.18)

The diﬀerentiation matrix 2.16 is the standard one in the literature of spectral methods [6, 25, 33]. The matrix actually used in the algorithm is a modiﬁed version of it. In fact, the CGL points are usually deﬁned from 1 to −1, but for control optimization problems it is more convenient to deﬁne them from −1 to 1 as in eqn. (2.10). In this case, the correct diﬀerentiation ˜ N is: matrix D ˜ N = −DN D (2.19) as described in [10]. For the second derivative calculation, the second-order diﬀerentiation ma2 trix D˜N is used: 2 2 ˜ 2 β. ¨ β= D (2.20) N tend − tini

During the approach phase, the thrust is reduced to gain manoeuvrability; the discontinuity in the thrust proﬁle can’t be accurately represented through a pseudospectral approximation. The simulation has been splitted in 2 domains to handle control discontinuities. • Simulation phase 1 (SYM1) It covers the Main Brake and the ﬁrst part of the Approach phase. The main thrust has a constant value. 0 ≤ t ≤ t1 , T1min ≤ T1 (t) = T1 ≤ T1max , −π/2 ≤ β1k ≤ π/2, k = 0, 1, . . . , N1β .

(2.21) (2.22) 13

Chapter 2 • Simulation phase 2 (SYM2) It covers the second part of the Approach phase till the Low Gate, with throttleable thrust. t1 ≤ t ≤ tf , T2min ≤ T2k ≤ T2max , k = 0, 1, . . . , N2T , −π/2 ≤ β2k ≤ π/2, k = 0, 1, . . . , N2β .

(2.23) (2.24)

The system states are imposed to be continuous between the two simulation domains: the resulting states at the end of SYM1 are used as initial condition of SYM2. The lander attitude cannot present discontinuities; this implies the additional control constraint: β1N1β = β2 0 .

(2.25)

After the High Gate, the trajectory must be compatible with navigation and Hazard Detection and Avoidance (HDA) constraints, in terms of visibility of the NLS and time for image processing. Then, after the hazard detection, the lander must have the necessary divert capabilities, in terms of altitude and speed. This can be obtained by imposing additional constraints on states and control variables at the end of SYM1: hmin ≤ h(t1 ) ≤ hmax , u min ≤ u(t1 ) ≤ umax , (2.26) v ≤ v(t ) ≤ v , min 1 max βmin ≤ β1N1β ≤ βmax .

Also the maximum allowed control torque has a crucial role in the determination of the trajectory in the proximity of the landing site. In particular, a lower torque corresponds to a more vertical approach towards the NLS. The resulting problem is: Find the control profiles: • T1 ; • β1k ,

k = 0, 1, . . . , N1β ;

• T2k ,

k = 0, 1, . . . , N2T ;

• β2k ,

k = 0, 1, . . . , N2β ;

and the maenuver times: t1 and tf that, applied to the system (2.1), with initial conditions (2.2), minimize the fuel consumption ffuel , according to control constraints (2.3), (2.5), (2.22), (2.24), (2.25), (2.7), (2.21), (2.23) and to state constraints (2.4), (2.26) and (2.6). 14

Nominal Landing Trajectory

2.2.2

Nominal Landing Simulation Results

R . The dynamic system (2.1) The problem has been implemented in Matlab has been integrated through a Runge-Kutta (4,5) method (ode45 ) and a nonlinear programming solver (SNOPT) has been used for the optimization. The physical constant values, and the assumptions on the lander architecture are summarized in Table 2.1. For a more detailed dissertation about the assumptions on the lander conﬁguration, see ref. [14, 15, 7] and the section 3.4 at page 34.

Table 2.1: Nominal landing simulation assumptions.

Physical Constants 4.9 × 1012 m3 s−2 1.738 × 106 m 9.806 65 m s−2

µ RM g0 Lander architecture

325 s 1000 kg m2 1000 N 2320 N 2500 N 3800 N

Isp Imax T1min T1max T2min T2max

0 −10

4000

Sym phase 1 Sym phase 2

3500

−30

3000

Thrust Reduction

Thrust [N]

Angle of Thrust [deg]

−20

−40 −50 −60

2500

2000

HG AG

−70

1500

−80 −90 0

Sym phase 1 Sym phase 2

100

200

300 Time [s]

400

(a) Angle of thrust

500

600

1000 0

100

200

300 Time [s]

400

500

600

(b) Thrust magnitude

Figure 2.3: Nominal landing path. Control variables profiles.

15

Chapter 2 In order to make the trajectory compatible with hazard detection requirements, a trade-oﬀ on constraints (2.26) imposed at the end of the SYM1 phase, and on the maximum allowed control torque (2.7) has been made. A time of 10 s needed by the Hazard Detection System to elaborate images is assumed [24]. A further margin of safety of 2 is applied, for a total of 20 s. The most suitable values resulted are summed up in Table 2.2. Table 2.2: Nominal landing simulation - Additional constraints values.

Constraint h(t1 ) u(t1 ) v(t1 ) β(t1 )

Min

Max

2000 −30 −35 −35

U.o.m. m m s−1 m s−1 deg

2500 30 0 35

Maximum control torque Mmax

Nm

30

The chosen orders of approximation are: N1β = 11, N2T = 9, N2β = 5. The resulted control proﬁles are shown in Figure 2.3, in terms of angle and magnitude of thrust. The thrust angle β is linked to the pitch angle θ by the relation θ = −β − π. During the Main Brake phase the optimal thrust value is the maximum available, while a bang-bang proﬁle appears with throttleable thrust. Corresponding states evolution is shown in Figures 2.4 and 2.5; the overall trajectory is depicted in Figure 2.7 at page 18. 20

2000

Horizontal speed [m/s]

Vertical speed [m/s]

0

−20

−40

−60

−80 0

100

200

300 Time [s]

400

(a) Vertical Speed

500

600

1500

1000

500

0 0

100

200

300 Time [s]

(b) Horizontal Speed

Figure 2.4: Nominal landing path. Speed trend.

16

400

500

600

Nominal Landing Trajectory 2

x 10

4

1500 1400

1.5 Mass [kg]

Altitude [m]

1300

1

1200 1100 1000

0.5

900

0 0

100

200

300 Time [s]

400

500

600

800 0

100

200

(a) Altitude

300 Time [s]

400

500

600

(b) Mass

Figure 2.5: Nominal landing path. Altitude and mass trend.

Figure 2.6 shows the trend of the view angle onto the NLS during the descent. The view angle is deﬁned as the angle between the vertical on the NLS and the spacecraft position, and it can vary from 0 (vertical of the LS) to 90◦ (LS at the horizon). Values higher than 90◦ mean that the landing site is below the horizon and it is not in sight at all. 100 NLS under Horizon

90

NLS View Angle [deg]

80 70 60

HG

50 40

AG

30 Thrust Reduction

20 10 0 0

Sym phase 1 Sym phase 2 100

200

300 Time [s]

400

500

600

Figure 2.6: View angle onto Nominal Landing Site.

The High Gate is deﬁned as the moment that NLS comes into the navigation sensors Field of View, and is assumed as the moment in which the view angle takes a value of 60◦ . The Approach Gate is intended as the instant in which the visibility onto the landing site area respects the operative require17

Chapter 2 ments of the Hazard Detection system, whose performances are determined by the trajectory inclination and by the distance from the target [16]. Then, is considered that at AG the NLS must be visible with an angle of view not greater than 45◦ , at an altitude not higher than 3000 m. In order to maximize the divert capability, the desired altitude at which the ﬁrst retargeting is possible is at least 2000 m. From Figure 2.6 is visible that an angle of 45◦ is reached 500.9 s after PDI; from Figure 2.5a it can be seen that the altitude is 3000 m at t = 511.2 s and 2000 m at t = 532.6 s. The time available for the Hazard Detection system is at least 21.4 s. 6

x 10 1.75

PDI

Main Brake

HG

1.7 Trajectory Moon Surface

1.65 −5

−4

−3

−2

−1

0 5

x 10 (a) Complete trajectory 6

1.676

x 10

AG

1.675 Thrust Reduction 1.674 1.673 1.672 1.671 −4.82

−4.8

NLS −4.78

−4.76

−4.74

−4.72

(b) Zoom on approach phase

Figure 2.7: Nominal landing path. Overall trajectory.

18

−4.7 5 x 10

Chapter

3

Retargeting Algorithm In this chapter the complete retargeting algorithm development is explained. First, all the reference systems and the equations of motion used are listed. Then, the trajectory modelling and the subsequent fuel optimization algorithm adopted are expounded. Performance and precision obtained are analysed through Monte Carlo method. Finally, it is showed that the solution found by the proposed algorithm is an approximation of the optimal one.

3.1

Landing Model

The retargeting problem, as part of an hazard detection and avoidance system, involves only the last part of the landing. Distances, for both downrange and altitude, are small (less than 3 km) compared to the planet’s radius and the assumption of a constant gravity ﬁeld with ﬂat ground is appropriate. Dealing with a lunar landing, no aerodynamic forces are considered. The eventual presence of atmosphere (especially with low density, as in the case of Mars) could be however negligible, due to the relative low velocity involved (on the order of 100 m s−1 at the beginning of the approach phase). In this case these forces, omitted in the trajectory design, can be treated as disturbances by the control system [2].

3.1.1

Reference systems

Diﬀerent reference systems are needed for the problem analysis. The relationship between them is shown in Figure 3.1. • Body-Fixed Indicated with (xB , yB , zB ), is centered on the lander center of mass, 19

Chapter 3 and aligned with the principal axes of inertia, as shown in Figure 3.2 • Flight Axes The ﬂight reference system (xF , yF , zF ) is centered on the lander center of mass, with the zF axis pointing towards the ground, and the xF axis lying in the plane of the nominal landing path, towards the NLS. The lander attitude is expressed as the rotation of the body-ﬁxed reference system from ﬂight axes. • Ground Axes The ground reference system (xG , yG , zG ) is centered onto the landing site. xG is zenith-pointed, while yG is aligned with xF . As the target landing site changes (consequently to a retargeting), the reference system moves with it without rotating.

Figure 3.1: Body-Fixed, Flight and Ground reference systems.

3.1.2

Equation of motions

In the trajectory design process, only the translational dynamics of the lander is considered. No thrust vectoring is considered, and the main thruster (or cluster of thrusters) is assumed to be tightly linked to the lander body. The dynamic of the center of mass, expressed in ground reference system, can be described by a set of 7 ordinary diﬀerential equations (ODE): • 3 are relative to the position vector rG = {rxG , ryG , rzG }T ; 20

Retargeting Algorithm

Figure 3.2: Body-Fixed reference system and euler angles.

• 3 for the speed vector vG = {vxG , vyG , vzG }T ; • 1 for the lander mass m. The acceleration of the center of mass is provided by the main thrust vector, written in ground axes, TG = {TxG , TyG , TzG }T , and by the gravity ﬁeld contribution gG = {−gM , 0, 0}T . Eqns. (3.1) show the resulting ODE system, where T = kTk is the thrust magnitude, Isp is the speciﬁc impulse of the main thrust, and g0 is the standard gravity acceleration: r˙ G = vG , G TG + gG , v˙ = (3.1) m T . ˙ =− m Isp g0

Since the main thruster is assumed to be rigidly connected to the spacecraft, the direction of the thrust vector is determined directly by the lander attitude. This can be expressed in Euler angles, ea = {φ, θ, ψ}T , in the 231 form. Figure 3.2 shows how angles are considered: φ is the Roll angle, θ is the Pitch angle and ψ is the Yaw angle. The 231 form is preferred to the more traditional 321, because allows to avoid the presence of singularities in the angles determination, between the boundaries imposed in the approach phase 21

Chapter 3 (see section 3.2). The eqn. (3.2) shows how the attitude director cosine matrix is obtained by the Euler angles. The lander attitude can be expressed also by a unit quaternion vector q = {q1 , q2 , q3 , q4 }T . In this case, the director cosine matrix can be written in the form (3.3): cψcθ sψ −cψsθ cφsψsθ + sφcθ , (3.2) ADC = −cφsψcθ + sφsθ cφcψ sφsψcθ + cφsθ −sφcψ −sφsψsθ + cφcθ ADC

2 q4 + q12 − q22 − q32 2(q1 q2 + q3 q4 ) 2(q1 q3 − q2 q4 ) 2(q2 q3 + q1 q4 ) . = 2(q1 q2 − q3 q4 ) q42 − q12 + q22 − q32 2 2(q1 q3 + q2 q4 ) 2(q2 q3 − q1 q4 ) q4 − q12 − q22 + q32

The rotation of the ﬂight reference system and deﬁned by the matrix: 0 0 AGF = 1 0 0 −1

(3.3)

from the ground axes is constant −1 0. 0

(3.4)

The thrust vector in Body axes is ﬁxed, varying only in magnitude: TB = T {−1, 0, 0}T . Its components in Ground axes can be obtained as: TG = ATGF ADC T TB

(3.5)

The system (3.1) can be rewritten, expressing the thrust vector with Euler angles. Into eqns. (3.6), vectors are expanded, and rG , vG , TG are indicated with r, v and T for a more clear notation: r˙x = vx , r˙y = vy , r˙z = vz , cos ψ sin θ v˙ x = −T − gM , m (3.6) cos ψ cos θ v˙ y = −T , m sin ψ v ˙ = T , z m T m . ˙ =− Isp g0 22

Retargeting Algorithm Using quaternions: r˙x = vx , r˙y = vy , r˙z = vz , 2 q q − q q 1 3 2 4 − gM , v˙ x = T m q22 + q32 − q12 − q42 v ˙ = T , y m 2 q q + q q 1 2 3 4 , v˙ = T z m T m . ˙ =− Isp g0

3.2

(3.7)

Trajectory Design

As shown by Equations (3.6), the system has 7 states variables (r, v and m); in this notation the Euler angles θ and ψ, together with the thrust magnitude T are seen as control variables. The roll angle φ doesn’t aﬀect the landing and it is always considered to be zero. As retargeting is requested, the maneuver starts at time t0 and ends at time tf . tToF = tf − t0 is deﬁned as time-offlight. All state and control variables with the subscript 0 are intended as their values at time t0 . The subscript f has the same meaning for the time tf . A correct landing trajectory has to satisfy a set of boundary constraints: • Initial position and velocity vectors (state constraints): r0 = {rx0 , ry0 , rz0 }T , v0 = {vx0 , vy0 , vz0 }T ; • Final desired position and velocity vectors (state constraints): rf = {rxf , ryf , rzf }T , vf = {vxf , vyf , vzf }T ; • Initial acceleration vector (control constraints): a0 = {ax0 , ay0 , az0 }T ; • Final acceleration vector (control constraints): af = {axf , ayf , azf }T ; 23

Chapter 3 At the end of the approach phase, a soft landing requires the spacecraft to be onto the vertical of the landing site, at 30 m altitude. The horizontal speed must be near zero, and a vertical speed of −1.5 m s−1 is required [15]. Then: rf = {30, 0, 0}T ,

(3.8) T

vf = {−1.5, 0, 0} .

(3.9)

Acceleration constraints are classiﬁed as control constraints. In fact, from Eqns. (3.6) can be seen that they are functions of control variables: cos ψ sin θ ax = −T − gM , m cos ψ cos θ (3.10) ay = −T , m az = T sin ψ . m

At the maneuver start, the lander attitude is known and ﬁxed; the acceleration depends only by the thrust magnitude at initial time (T0 ). On the other hand, at maneuver’s end the lander is desired to have a vertical attitude. The corresponding values of Euler angles are: θ(tf ) = −π/2, ψ(tf ) = 0. With these values, the acceleration vector at ﬁnal time is reduced to af = {Tf /m − gm , 0, 0}T . In this case, in order to minimize the number of free variables, the value of Tf is left free. Then, the acceleration vector at tf gives only 2 boundary constraint on axes y and z: ( af y = 0, (3.11) af z = 0. For a total of 17 boundary constrains: 6 for initial known boundary conditions on r0 and v0 ; 3 on a0 expressed by Eqns. (3.10), 6 for ﬁnal boundaries on position rf (3.8) ans speed vf (3.9), and 2 for the ﬁnal acceleration components aﬀected by the vertical landing constraint (3.11).

3.2.1

Polynomial Formulation

The trajectory can be expressed, separately for each axis, in a polynomial form, of minimum degree necessary to satisfy boundary constraints. 24

Retargeting Algorithm There are 6 boundary constraints for each horizontal axis (y and z): initial and ﬁnal acceleration, velocity and position. The position versus time proﬁle needed is of order 5. For a more clear notation, is convenient to start from the acceleration proﬁle represented by a cubic equation:

ay (t) = a0y + C1y t + C2y t2 + C3y t3 , az (t) = a0z + C1z t + C2z t2 + C3z t3 .

By integrating two times, and immediately imposing initial conditions at t0 :

1 1 1 vy (t) = v0y + a0y t + C1y t2 + C2y t3 + C3y t4 , 2 3 4 1 1 1 vz (t) = v0z + a0z t + C1z t2 + C2z t3 + C3z t4 , 2 3 4 1 1 1 1 ry (t) = r0y + v0y t + a0y t2 + C1y t3 + C2y t4 + C3y t5 , 2 6 12 20 1 1 1 1 rz (t) = r0z + v0z t + a0z t2 + C1z t3 + C2z t4 + C3z t5 . 2 6 12 20

There are only 5 boundary constraints for the x component. This leads to a quadratic acceleration proﬁle:

ax (t) = a0x + C1x t + C2x t2 , 1 1 vx (t) = v0x + a0x t + C1x t2 + C2x t3 , 2 3 1 1 1 rx (t) = r0x + v0x t + a0x t2 + C1x t3 + C2x t4 . 2 6 12

Now, by imposing the ﬁnal boundary conditions, constant terms can be ob25

Chapter 3 tained: rf x r0x vf x v0x a0 x , − 24 3 − 6 2 − 18 2 − 6 3 tf tf tf tf tf r0x v0x vf x a0x rf x = −36 4 + 36 4 + 24 3 + 12 3 + 6 2 , tf tf tf tf tf rf y a0y r0y vf y v0y af y = 60 3 − 60 3 − 24 2 − 36 2 + 3 −9 , tf tf tf tf tf tf r0y vf y v0y af y a0y rf y = −180 4 + 180 4 + 84 3 + 96 3 − 12 2 + 18 2 , tf tf tf ff tf tf rf y r0y vf y v0y af y a0y = 120 5 − 120 5 − 60 4 − 60 4 + 10 3 − 10 3 , tf tf tf tf tf tf a0z r0z vf z v0z af z rf z −9 , = 60 3 − 60 3 − 24 2 − 36 2 + 3 tf tf tf tf tf tf rf z r0z vf z v0z af z a0z = −180 4 + 180 4 + 84 3 + 96 3 − 12 2 + 18 2 , tf tf tf ff tf tf r0z vf z v0z af z a0z rf z = 120 5 − 120 5 − 60 4 − 60 4 + 10 3 − 10 3 . tf tf tf tf tf tf

C1x = 24 C2x C1y C2y C3y C1z C2z C3z

The trajectory is completely deﬁned, function of tf and T0 . While tf directly inﬂuences the constants values, T0 does the same indirectly, determining the initial acceleration components through the eqn. (3.10). Then, a complete guidance proﬁle can be extrapolated. The thrust-to-mass ratio vector P(t) is simply obtained from the acceleration proﬁle. From eqn. (3.1): ax (t) + gM T(t) ay (t) . P(t) = = (3.12) m az (t) Then, from eqn. (3.6) the mass versus time trend is determined by a 1st order linear ordinary diﬀerential equation. If P (t) = kP(t)k: P (t) m(t), Isp g0 Z m(t) = m0 exp −

(3.13)

m(t) ˙ =−

tf

t0

! P (t) dt . Isp g0

(3.14)

The analytical calculation of the integral exponent is complex, but it can be easily obtained trough numerical integration, using Chebyshev pseudospectral methods, such as the Clenshaw-Curtis quadrature [33], through 26

Retargeting Algorithm which the integral is discretized to a ﬁnite sum: Z tend Z 1 tend − tini 1 1 P (t)dt = − − P (τ )dτ Isp g0 tini Isp g0 2 −1 N 1 tend − tini X wCk Pk , =− Isp g0 2

(3.15)

k=0

where Pk are the values of the thrust-to-mass ratio at CGL points (see chapter 2.2), and wC are the weights of the Clenshaw-Curtis quadrature scheme (see Ref. [33], chapter 12). The system mass can be evaluated at every time instant t, taking tini = t0 , and tend = t. The choice of the approximation order N aﬀects both the precision and the computational speed of the algorithm. The topic is discussed in chapter 3.4. Through the mass versus time trend it is possible to obtain the thrust vector T and its magnitude T from the thrust-to-mass ratio: T = P/m.

(3.16)

ˆ can be exploited to obtain the lander’s attitude in The thrust unit vector n Euler angles: − cos ψ sin θ T ˆ= n = − cos ψ cos θ , (3.17) T sin ψ n ˆx , −π ≤ θ ≤ 0, θ = arctan n ˆ y n ˆz (3.18) , −π/2 ≤ ψ ≤ π/2, ψ = arctan p 2 2 n ˆx + n ˆy φ = 0.

A discrete guidance proﬁle, in terms of thrust magnitude and Euler angles at CGL points, is obtained, function of T0 and tToF . Guidance can be reconstructed at every time instant, with the polynomial approximation 2.14.

3.2.2

Additional Constraints

The problem is so reduced to ﬁnd the values of T0 and tToF , according to any additional constraint not implicitly satisﬁed by the polynomial formulation, in order to minimize the fuel consumption. 27

Chapter 3 Assuming x = {tToF , T0 }T , the cost function is f (x) = m0 − mf , and the problem can be expressed in the form: ( xL ≤ x ≤ xU min f (x) such that . (3.19) x cL ≤ c(x) ≤ cU The optimization variables x are not allowed to take any value, but they have a ﬁnite domain with lower bound xL and upper bound xU . These are called box constraints. The elements of c(x) are generally non-linear functions of the optimization variables, also bounded between lower and upper limits cL and cU . These constraints need to be satisﬁed during all the landing maneuver, so they are called path constraints. Finally, the polynomial formulation does not explicitly consider boundary constraint on mass. This implies the additional constraint: (3.20)

mdry ≤ mf ≤ m0 .

This constraint is implicitly satisﬁed at its upper bound, but not at its lower one. It it equivalent to impose an upper limit to the cost function. Box Constraints The thrust magnitude is bounded to the thrust available on-board: 0 < Tmin ≤ T0 ≤ Tmax .

(3.21)

The time-of-ﬂight must be greater than zero, while its theoretical upper limit is determined by the maximum amount of fuel on board mfuel . 0 ≤ tToF ≤ tmax = mfuel

Isp g0 . Tmin

(3.22)

Path Constraints • Thrust Bounds As the box constraint on T0 , during all the landing phase the required thrust magnitude cannot exceed the limit imposed by the actual engine on board: Tmin ≤ T (t) ≤ Tmax . (3.23) • Control Torques Euler angles rate of change is subject to the actual control torques available by the ACS. The extrapolation of the exact torques from angles is not immediate, due to coupled terms in the attitude dynamics. 28

Retargeting Algorithm The objective is to characterize such a rotational rate constraint without coupling the problem to the rotational dynamics, in order to save computation time. Torques are approximated by the decoupled term due to the angular acceleration. This corresponds to the reality in case of small angles and low angular speed, a condition not at all veriﬁed during the approach phase, and it can be used only for an estimate: ¨ MCy = Iy θ, ¨ MCz = Iz ψ, where Iy and Iz are the principal moments of inertia on y and z axes (body axes). The actually considered constraints has the form: ¨ ≤ ρMCmax , −ρMCmax ≤ Imax θ(t) ¨ ≤ ρMCmax . −ρMCmax ≤ Imax ψ(t)

(3.24) (3.25)

As the fuel is expended, moments of inertia decrease with the lander mass. Imax is the maximum moment of inertia at initial time t0 . This allows to avoid the on-board calculation of inertia properties, and gives a margin of safety in the torques calculation. An additional margin of safety 0 < ρ < 1 can be applied. • Glide-slope constraint In a feasible landing path altitude is always greater than zero. This constraint can be improved considering a glide-slope constraint. In this case the lander is required to remain in a cone deﬁned by the minimum slope angle γmax , as showed in Figure 3.3. This constraint has a dual purpose: it assures that the the lander does not penetrate the ground, even in presence of bulky terrain features near the landing site; at the same time it limits the angle of view onto the target. In fact, the performances of vision-based navigation systems depend on inclination between the trajectory and the ground [11, 27]. Every state constraint on position and velocity can be expressed in the general form [2]: kSj w(t) + dj k + cTj w(t) + qj ≤ 0,

(3.26)

where Sj ∈ Rnj ×6 , dj ∈ Rnj , with nj ≤ 6, cj ∈ R6 , qj ∈ R and w ∈ R6 is the part of the state vector corresponding to position and velocity deﬁned by: r(t) . (3.27) w(t) = v(t) 29

Chapter 3 Position Trajectory

γmax Glide-slope constraint

Target LS

Figure 3.3: Glide-slope constraint. This constraint requires to the spacecraft to remain in a cone defined by the minimum slope angle γmax . The apex of the cone coincides with the target landing site.

For the glide-slope constraint, in order to have 0 ≤ γ(t) ≤ γmax ≤ π/2: 0 1 0 0 0 0 S= , 0 0 1 0 0 0 cT = − tan γmax 0 0 0 0 0 , d = q = ∅.

The constraint take the form: − ∞ ≤ kSw(t)k + cT w(t) ≤ 0.

(3.28)

Path constraints need to be satisﬁed at every time instant during the landing. The algorithm evaluates them discretely, at CGL points. In the case of the control torques constraint, a second derivative of the Euler angles is necessary. It is obtained by the use of the Chebyshev diﬀerentiation matrix with the same method expounded in chapter 2.2.

3.3

Optimization Algorithm

The optimization problem (3.19) can be solved through many diﬀerent algorithms. In this case, fast computation must be privileged, in perspective of a real-time implementation for on-board hardware. In this context, direct optimization methods are attractive, because they don’t require any derivation of necessary conditions, treating the cost function as a “black-box” [10, 1]. 30

Retargeting Algorithm Preliminarily, a modiﬁed version of compass search, enhanced in order to handle also non-linear constraints, is proposed. The compass search method, also known as coordinate search, belongs to a broad class of derivative-free algorithms known in literature by the name of Generating Set Search (GSS) methods [20]. A GSS method samples the objective function f , along a suitable set of search directions, around the current iterate xk , at a given distance depending on a stepsize parameter ∆k . In a classical coordinate search method, the search directions consist of the coordinate axis, and the new iterate xk+1 is sought in the set of trial points P + {xk ± ∆k e1 , xk ± ∆k e2 , . . . , xk ± ∆k en }. If P contains a point satisfying a suitable decrease condition, such point is accepted as the new iterate. Otherwise, the stepsize is reduced. The algorithm stops when the stepsize ∆k falls below a prescribed threshold τ∆ . First, the optimization variables are normalized, in order to give them the same relative weight in the optimization: ˜= x

x − xL xU − xL

⇔

˜ (xU − xL ) + xL . x=x

(3.29)

This lead to a modiﬁcation of the box constraints: ˜ ≤ 1. 0≤x

(3.30)

Then a feasibility function F (˜ x) is created, deﬁned as: NC X 1 F (˜ x) = max(0, c˜j ), wF j j=0

(3.31)

where c˜j are the components of a generalized constraints vector ˜c(˜ x): cL − c(˜ x) c(˜ x) − cU , (3.32) c˜(˜ x) = 0−x ˜ ˜−1 x

and wF is a vector of weights, used to normalize diﬀerent constraints that can have diﬀerent orders of magnitude: cU − cL cU − cL (3.33) wF (x) = xU − xL . xU − xL 31

Chapter 3 The glide-slope lower bound (3.28), and consequently the corresponding weight, is inﬁnite. In order to avoid an improper constraint evaluation, this weight is setted to r0x , that have the correct order of magnitude. ˜ leads to a feasible trajecIn this way, if a set of optimization variables x tory, the corresponding feasibility function has value 0. On the contrary, in case of infeasibility F (˜ x) > 0. The optimization algorithm, called GuidALG, operates in two phases: in ˜ 0 , it searches for any feasible set the ﬁrst one, starting from an initial point x of optimization variables; in the second one, it seeks in the feasible domain for an optimal solution. Feasibility Step ˜ 0 at The feasibility step ﬁrstly checks the feasibility of the starting point x the center of the optimization variables domain. If feasible, pass directly to the optimality step; if not, it performs an unconstrained compass search on the function F (˜ x). The search is stopped when a feasible point is found (F (˜ x) = 0), or when the iteration limit is reached. In this case, the problem is classiﬁed as infeasible. The procedure is described in the ﬁrst part of algorithm 1. Optimality Step The optimality step inherits the starting point and the mesh size ∆ from the feasibility step, if successful. Then, it performs an unconstrained compass search on the modiﬁed cost function Φ(˜ x), deﬁned as: Φ(˜ x) = f (˜ x) + η sgn F (˜ x) , (3.34)

where f (˜ x) is the original cost function of the problem (3.19), and η = 100 10 , a number that is certainly greater than the maximum value that cost function can have. The search is stopped when the mesh size decrease under a minimum value τ∆ , or when the maximum iteration limit is reached, as described in the second part of algorithm 1.

32

Retargeting Algorithm

Algorithm 1 [x, fea, opt] ← GuidALG(x0 ∈ R2 , ∆0 > 0, τ∆ , itF max , itOmax ) 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24:

25: 26: 27: 28: 29: 30:

k ← 0; fea ← TRUE; ⊲ Initial point feasibility if F (˜ x0 ) > 0 then ⊲ Infeasible x˜0 : Feasibility Step Start while k < itF max do x⋆ ← argmint F (t) , t ∈ Pk + {˜ xk ± ∆k ei : i = 1, 2}; if F (x⋆ ) = 0 (feasible point) then x˜k+1 ← x⋆ ; k ← k + 1; break; ⊲ Feasibility - go to next step else if F (x⋆ ) < F (xk ) then x˜k+1 ← x⋆ ; ∆k+1 ← ∆k ; k ← k + 1; else ∆k+1 ← ∆2k ; k ← k + 1; end if end while if k = itF max then fea ← FALSE; return; ⊲ Infeasible - Return end if end if ⊲ Feasibility Step End j ← k; opt ← TRUE; ⊲ Optimality Step Start while j < itOmax + k do xj ± ∆j ei : i = 1, 2}; x⋆ ← argmint Φ(t) , t ∈ Pj + {˜ if Φ(x⋆ ) < Φ(˜ xj ) then x˜j+1 ← x⋆ ; ∆j+1 ← ∆j ; j ← j + 1; else if ∆j ≤ τ∆ (minimum mesh size) then j ← j + 1 return; ⊲ Optimality reached else (mesh local optimizer ) ∆ ∆j+1 ← 2j ; j ← j + 1; end if end while if j = itOmax + k then opt ← FALSE; return; ⊲ Max iteration reached end if ⊲ Optimality Step End

33

Chapter 3

3.4

Algorithm Performances

Algorithm performances have been evaluated through Monte Carlo simulations, in order to estimate computation time, divert capabilities, and optimality of the solution found. The earliest retargeting altitude considered is 2000 m. R code, and tested The retargeting algorithm has been written in Matlab R

TM on a Intel Core i7-2630QM CPU at 2 GHz of frequency. As practical example, the ESA Lunar Lander mission is taken as reference. The Lunar lander proposed propulsion system is made by ﬁve 500 N EAM main engines, only 2 of them active during the approach phase, and by six 220 N ATV thrusters, pulse-modulated in order to obtain throttleable thrust between 1000 N and 2320 N [14, 15]. In the real case, the mixed use of diﬀerent engines implies a varying speciﬁc impulse, considered as constant in simulations. Also canted out by ϑP

ϑP

Roll Thrusters Pitch Thrusters Yaw Thrusters

z x Figure 3.4: Considered RCS thrusters scheme.

A Reaction Control System made by 4 clusters of Astrium ST 22 N bipropellant thrusters is assumed [15]. Figure 3.4 shows the considered RCS scheme. Supposing the cant angle θP = 15◦ and lander dimensions showed in ref. [7] (see also Figure 3.5), this lead to a control torque of 54.4 N m on pitch and yaw axes; a value of 50 N m has been adopted in simulations. The maximum moment of inertia Imax has been calculated considering the right cylinder of the lander central body, as in Figure 3.5, assuming the mass at Powered Descent Initiation (PDI, see section 2.1) at page 7). The resulting value is certainly greater than the mass at ﬁrst retargeting, giving, along with the reduced maximum control torque considered, a margin of 34

Retargeting Algorithm

Figure 3.5: ESA Lunar Lander dimensions (adopted from ref. [7]).

safety of about 2 on the control torques constraint evaluation. The lander mass at 2000 m has been extrapolated from simulations of nominal trajectory (see section 2.2 at page 9). All the assumptions about the lander architecture are summarized in Table 3.1. A glide-slope constraint of 70◦ has been considered. The standard gravity acceleration at sea level is g0 = 9.806 65 N m−2 and the gravity acceleration on Moon’s surface is gM = 1.624 681 N m−2 . Table 3.1: Lander architecture assumptions.

Feature

Value

Mass @ PDI Mass @ 2000 m altitudea Dry mass Isp Imax Tmin Tmax MCmax

1500 kg 865 kg 790 kg 325 s 1000 kg m2 1000 N 2320 N 50 N m

a

3.4.1

From simulation, see chapter 2.2

Order of approximation

The polynomial formulation of the landing trajectory produces inherently exact states at the maneuver end. On the other hand, the thrust proﬁle 35

Chapter 3 is computed through numerical methods, due to integral term of the mass calculation. The choice of the order of approximation N aﬀects both precision and computational speed of the algorithm: higher degree improves the precision of the path evaluation, but slows the computation (increasing the calculation time of single iterations). The computational cost of the algorithm has been evaluated by Monte Carlo simulations. Infeasible points are not considered in this analysis, since the calculation time of them depends essentially by the maximum number of feasibility iteration itF max allowed, a parameter which requires a tradeoﬀ itself. A set of 100 000 feasible points for each value of N has been used. The retargeting is assumed to be ordered from the nominal landing path, at a random altitude between 100 m and 2000 m, with a random ordered diversion between -2000 m and +2000 m in both the horizontal directions. Results are shown in Figure 3.6. 80 mean 3σ

Computation Time [ms]

70

60

50

40

30

20 10

15

20 N

25

30

Figure 3.6: Mean computation time varying order of approximation N

In order to estimate the algorithm precision, the control proﬁle calculated by the algorithm, reconstructed with the eqn. (2.14), has been applied to the system (3.6) and integrated trough a traditional Runge-Kutta (4,5) method. The result is compared to the desired one to obtain the error on ﬁnal position and speed. For each value of N, a set of 10 000 feasible points has been considered. In order to obtain the worst case, the initial altitude is constant at 2000 m, in order to maximize dispersions due to the approximation of the guide proﬁle 36

Retargeting Algorithm reconstruction. The commanded diversion is random between -2000 m and +2000 m in both the horizontal directions. Results are shown in Figure 3.7 for the position error, in Figure 3.8 for the velocity. 0

0

10 Vertical Position Error [m]

Horizontal Position Error [m]

10

−1

10

mean 3σ

−2

10

−3

10

−4

10

10

−1

10

mean 3σ

−2

10

−3

10

−4

15

20 N

25

10

30

10

15

(a) Horizontal.

20 N

25

30

(b) Altitude.

Figure 3.7: Position precision varying order of approximation N.

−2

−2

10 Vertical Speed Error [m/s]

Horizontal Speed Error [m/s]

10

−3

10

mean 3σ

−4

10

−5

10

−6

10

10

−3

10

mean 3σ

−4

10

−5

10

−6

15

20 N

25

(a) Horizontal.

30

10

10

15

20 N

25

30

(b) Altitude.

Figure 3.8: Speed precision varying order of approximation N.

It can be seen that the computation time is very stable, with very low dispersion. The polynomial approximation imposes on the trajectory a smooth proﬁle, on which the pseudospectral approximation is very eﬀective. The error is small, and it can considered negligible from N = 20 onwards. For all subsequent simulation, N is set to 20.

3.4.2

Divert Capability

Divert capability has been estimated from the nominal landing trajectory, at diﬀerent heights. Five values have been taken into considerations: 2000, 37

Chapter 3 1500, 1000, 500 and 100 m; for each value, 100 000 random target landing sites has been tested. As one might expect, the feasible diversion decreases with the altitude, as showed in Figure 3.9. In some conditions GuidALG, while converging towards an optimal solution, stops the optimization reaching the maximum allowed iteration number.

(a) r0x 2000 m; ±2000 m area.

(b) r0x 1500 m; ±2000 m area.

(c) r0x 1000 m; ±2000 m area.

(d) r0x 500 m; ±600 m area.

Figure 3.9: Feasibility performances. Divert capability from nominal landing path. Green: feasible; blue: infeasible; red: feasible, optimality search ended for maximum iteration limit reached; N = 20, itF max = 50, itOmax = 70. Graphs not at the same scale (Follows in Figure 3.10).

Analysing computational performances (see Figure 3.11a) it can be seen that the computation time out of the feasible area is practically constant, and depending, besides the approximation order N, from the maximum number of feasibility iterations allowed itF max (as showed in Figure 3.11b). 38

Retargeting Algorithm

Figure 3.10: Feasibility performances (follows from Figure 3.9). 100 m; ±100 m area.

Altitude r0x

Figure 3.12 shows the fuel consumption (that is the cost function of the optimum problem) from the altitude values of 2000 m and 100 m. Discontinuities and not continuous zones of feasibility indicate that the compass search algorithm in some cases can converge to suboptimal solution (see section 3.4.3).

(a) Computation time.

(b) Iteration number.

Figure 3.11: Computational performances. r0x 1000 m; ±2000 m area; N = 20, itF max = 50, itOmax = 70.

3.4.3

Optimality

The optimal solution of a landing problem has a well-known bang-bang thrust proﬁle [32, 13]. Three diﬀerent solutions of the problem have been compared, 39

Chapter 3

(a) r0x 2000 m; ±2000 m area.

(b) r0x 100 m; ±100 m area.

Figure 3.12: Fuel consumption performances.

in order to estimate the optimality achieved with the proposed method. • The solution computed by GuidALG (see algorithm 1 at page 33); • The proposed polynomial formulation optimized with a non-linear programming solver (SNOPT); • The solution of a direct pseudospectral collocation algorithm applied to the original problem (3.6), capable to handle discontinuous control proﬁles. Four representative cases are presented: • CASE 1: A recalculation of the trajectory pointing to the Nominal Landing Site, from 2000 m altitude (Figure 3.13); • CASE 2: An intermediate retargeting from 500 m of altitude, with a diversion of {0, +100, +100}T meters from the NLS (Figure 3.14); • CASE 3: A retargeting from 1000 m altitude, toward a point marked as feasible, but non optimal due to iteration limit reaching by GuidALG, at {0, −84.45, −277.46}T from NLS (Figure 3.16); • CASE 4: A diversion of {0, −1600, 0}T m from altitude 2000 m, a case of suboptimal solution by GuidALG. It is easy to see that the optimal solution of the polynomial formulation is an approximation of the real optimal trajectory, with consistent advantages of computational cost. In some cases, e.g. CASE 4, GuidALG converge towards a not optimal, but still feasible solution. 40

Retargeting Algorithm

2400

−40 θ [deg]

2200

−80 −100 0

1800 1600 1400

Direct Coll. Poly+SNOPT GuidALG

1200 1000 0

20

40

60

80

20

40 Time [s]

60

80

1

20

40 Time [s]

60

ψ [deg]

Thrust [N]

2000

−60

0

−1 0

80

(a) Thrust.

(b) Euler Angles.

860

830

855

829

850

Mass [kg]

Mass [kg]

865

845 840

827 826

835

825

830 825 0

828

20

40 Time [s]

(c) Mass.

60

80

824 60

62

64

66 Time [s]

68

70

(d) Mass - zoom.

Figure 3.13: Optimality comparison, CASE 1. Altitude 2000 m, TLS {0, 0, 0}T m. For this case, the compass search solution and the NLP solver solution coincide.

41

Chapter 3

2400

−60 θ [deg]

2200

−80

−100 0

1800 1600 1400

Direct Coll. Poly+SNOPT GuidALG

1200 1000 0

10

20

30

40

10

20 Time [s]

30

40

31

32

20

10

20 Time [s]

ψ [deg]

Thrust [N]

2000

30

0

−20 0

40

(a) Thrust.

(b) Euler Angles.

850

825.5

845 825 Mass [kg]

Mass [kg]

840 835

824.5

830 824

825 820 0

10

20 Time [s]

(c) Mass.

30

40

823.5 28

29

30 Time [s]

(d) Mass - zoom.

Figure 3.14: Optimality comparison, CASE 2. Altitude 500 m, TLS {0, 100, 100}T m.

42

Retargeting Algorithm

2400

−60 θ [deg]

2200

−80

−100 0

1800 1600 Direct Coll. Poly+SNOPT GuidALG

1200 1000 0

10

20

30

40

50

10

20 30 Time [s]

40

50

20

1400

10

20 30 Time [s]

ψ [deg]

Thrust [N]

2000

40

0 −20 −40 0

50

(a) Thrust.

(b) Euler Angles.

855 850

825.5

Mass [kg]

Mass [kg]

845 840 835

825

824.5

830

824

825 820 0

10

20 30 Time [s]

(c) Mass.

40

50

41

42

43

44

Time [s]

(d) Mass - zoom.

Figure 3.15: Optimality comparison, CASE 3. Altitude 1000 m, TLS T {0, −84.45, −277.46} m. In this case, GuidALG stops for maximum iteration limit reaching, while converging towards an optimal solution.

43

Chapter 3

2400

0 θ [deg]

2200

−100 −150 0

1800 1600 Direct Coll. Poly+SNOPT GuidALG

1200 1000 0

40

60

80

100

20

40 60 Time [s]

80

20

40 60 Time [s]

80

100

0

−1 0

100

(a) Thrust.

(b) Euler Angles.

870

826

860

824

850

822

Mass [kg]

Mass [kg]

20

1

1400

ψ [deg]

Thrust [N]

2000

−50

840

820 818

830

816 820 810 0

814 20

40 60 Time [s]

(c) Mass.

80

100

65

70

75 Time [s]

80

85

(d) Mass - zoom.

Figure 3.16: Optimality comparison, CASE 4. Altitude 2000 m, TLS {0, −1600, 0}T m. In this case, GuidALG converge toward a local minimum, achieving a suboptimal solution.

44

Chapter

4

Retargeting Simulation The guidance algorithm developed in chapter 3 has been tested through a simulation of the landing, from the ﬁrst retargeting possibility at 2000 m R altitude to the beginning of the TD phase. The model is written in Simulink environment. Actuators

Control

Guidance

Spacecraft Dynamics

Environment

Sensors

Navigation

HDA Manager

Figure 4.1: Complete landing system logical scheme.

Figure 4.1 illustrates the logical structure of a complete landing system equipped with an Hazard Detection and Avoidance system. On-board actuators and external environment inﬂuence the spacecraft dynamics with forces and torques. Information from inner dynamics and environment are transduced by sensors, and interpreted by the navigation block to obtain an estimate of the system states. The hazard detection manager cross-checks data from navigation and sensors to evaluate the safety of the landing site area, and commands a landing site change if necessary. With informations about the landing site and the current states, the guidance system computes 45

Chapter 4 the target attitude and thrust magnitude. Finally, these values are used by the control block to calculate the necessary commands for actuators.

4.1

Modeling criteria

In this work, not all the components of the system are modelled in detail: Figure 4.2 shows the simpliﬁed model adopted. The external environment aﬀects the lander dynamics through the gravitational ﬁeld and disturbance forces and torques. But external disturbance eﬀects have a very low order of magnitude, and they have a signiﬁcant impact only over relatively long times. Approach and Terminal Descent phases are characterized by a relatively short duration, in the order of magnitude of 100 s, and external disturbance eﬀects can be considered as negligible. The only environmental inﬂuence considered is the main gravitational ﬁeld, that is included in the lander dynamics. Sensors and navigation system are replaced by a navigation error model, that emulates errors in the states reconstruction process. A more detailed dissertation on errors generation is expounded in section 4.1.2. The model does not take into account the actual internal architecture of the Hazard Detection system. Predetermined diversions are commanded in order to evaluate guidance performances. The update rate considered for navigation, guidance and attitude controller during the Powered Descent phase is 20 Hz. Actuators

Control

Imposed Diversion

Spacecraft Dynamics Navigation Errors Model

Guidance

Figure 4.2: Simplified landing system.

4.1.1

Lander Dynamics

The lander translational dynamics is the same adopted in section 3.1.2, completed with rotational dynamics. The spacecraft is still assumed to be a rigid 46

Retargeting Simulation body. The dynamical system considered is: r˙ G = vG , TG G + gG , v ˙ = m 1 q˙ = Ωq, 2 ˙ B ), ω˙ B = I−1 (MBC + MBD − ω B × Iω B − Iω T ˙ =− , m Isp g0

(4.1)

where rG , vG , TG , gG , m, Isp and g0 have the same meaning of Eqns. 3.1. Superscripts G and B indicate a value expressed, respectively, in the Ground and in the Body-ﬁxed reference system, as deﬁned in section 3.1.1. The vector ω B = {ωxB , ωyB , ωxB }T contains the angular velocities. The vector of quaternions q represents the lander attitude ad deﬁned in section 3.1.2, while the matrix Ω is deﬁned as: 0 ωz −ωy ωx −ωz 0 ωx ωy . (4.2) Ω= ωy −ωx 0 ωz −ωx −ωy −ωz 0 MBC and MBD are, respectively, the control and the disturbance torques, in Body-ﬁxed frame, while I and I˙ are the matrix of inertia and its time derivative. For a more clear notation, superscript G and B are suppressed from this point forward. Lander Inertial Properties The lander inertial properties are obtained through a simple parametrization similar to the one adopted in [17]. In order to capture the main features with a simple model, basic geometric shapes with constant densities are used to approximate the various components. Since these shapes are all axially symmetric, the corresponding Centers of Mass (CoMs) all lie on the vehicle center-line. The locations of these CoMs are properly weighted by their respective masses to compute the location of the total CoM. Again, since the shapes are axially symmetric, they do not have products of inertia, only moments of inertia. These inertias are combined to compute the total inertia about the total CoM. Figure 4.3 47

Chapter 4

Figure 4.3: Parametric lander inertial model.

shows the division adopted. The mass allocation is summarized in Table 4.1. The ESA Lunar Lander, taken as reference, is assumed to have a mass of 2500 kg at launch [9], and a dry mass mdry = 790 kg is considered (see section 3.4 at page 34). Table 4.1: Components of lander inertial model.

Component

Shape

Top Structure Main Structure Fuel Bottom Structure

Circular Circular Circular Circular

Mass [kg] cylinder annulus cylinder cylinder

100 500 1710÷0 190

Combining the elementary properties of simple shapes, principal moments of inertia, polynomial functions of the lander mass, are obtained. Indicating through subscripts x, y and z the values relative to the body-ﬁxed axes as expounded in section 3.1.1: Ix = A0 + A1 m, Iy/z = B0 m−2 + B1 m−1 + B2 + B3 m + B4 m2 + B5 m3 .

(4.3) (4.4)

The time derivative is easily obtained at every time instant by: dI(m) dI(m) dm dI(m) T (t) . = =− dt dm dt dm Isp g0 48

(4.5)

Retargeting Simulation The resulting matrix of inertia at PDI (mass 1500 kg) is (in kg m2 ): 1204.7 0 0 I(1500) = 0 1070.0 0 , 0 0 1070.0

that becomes, when the lander reaches the dry mass mdry = 790 kg: 877.6 0 0 717.9 0 . I(790) = 0 0 0 717.9

(4.6)

(4.7)

Disturbance Torques Due to the relative small duration of the landing phase, external disturbance torques are not modelled. The only disturbance torque considered is due to the thrust misalignment during the descent. The point of application of the thrust vector T is assumed not aligned with the lander CoM, but translated along the body axes yB and zB . This shift is considered partly constant and partly random, due to the pulsed mode of the ATV thrusters assumed for throttleability (see section 3.4). The thrust point of application b in Body-ﬁxed axes can be expressed as: 0 0 b = by + ρb cos α , (4.8) ρ sin α bz b

where by and bz are the constant components of the shift, α is a uniform random angle between 0 and 2π, and ρb is a random distance of value 0 ± brand (1σ). The values for by , bz and brand are determined by some considerations on past missions and launcher requirements: • The Soyuz CSG launcher (a possible candidate for a lunar mission) requires that the CoM of the spacecraft stays within a distance d ≤ 15 mm from the launcher longitudinal axis [30]; • The Apollo 11 Lunar Module, a lander with a mass of one order of magnitude higher than the case under exam, had a misalignment between CoM and longitudinal axis equal to 22 mm [17]. The actual values assumed are: by = 5 × 10−3 m,

bz = −5 × 10−3 m, brand = 1 × 10−2 m. 49

Chapter 4 Since the thrust vector in Body-ﬁxed axes is T = T {−1, 0, 0}T , the resulting disturbance torque vector MD is: 0 0 0 MD = T 0 0 −1 b. (4.9) 0 1 0

4.1.2

Navigation Errors Model

Pinpoint landing requires high precision in the estimate of attitude, position and speed. During all the orbital phases, from the injection on the Lunar Transfer Orbit until PDI, the lander attitude is assumed to be estimated by one or more star trackers. During the Powered Descent phase greater angular rates and orientation changes can bring the lander out of star trackers operating range, and an Inertial Measurement Unit (IMU) is needed to keep the attitude under control. A visual-based navigation system is considered for the determination of position and speed. From the LLO and during the coasting phase, the absolute position is estimated through a comparison between the visible craters and those recorded in an on-board database [3, 29]. After the PDI, as the lander gets closer to the ground, the system passes to a relative navigation mode, in which the position and the speed are determined relatively to the ground, tracking landmarks during the descent. Attitude IMU directly measures rotational speeds through gyroscopes. The adopted error model is a simpliﬁed version of the one presented in [31] and in [34]. The measured angular velocities vector ω M is: 1+S M M ωx B nx M 1 + S M ω B ωM = + ny , (4.10) y + M M 1+S ωz B nz where S is the scale factor error, M the misalignment error, B the bias error ; nx , ny and nz are independent sensor noises, modelled as zero mean white noise. The attitude is obtained by integrating the angular speed, according to the relation: 1 q˙ = Ωq. 2 50

Retargeting Simulation This integration is assumed to be made on-board by the navigation system. A discrete time integrator, with update period 20 Hz is used. Is assumed that star trackers are disabled at PDI, and after this point the attitude determination relies only on gyroscopes propagation. Since that, the considered initial condition is the attitude assumed by the lander at PDI (q0 = {0, 0, 0, 1}T ), aﬀected by the star tracker error. The retargeting simulation starts at an altitude of 2000 m; before that, the attitude error is integrated from PDI along the nominal trajectory in order to obtain the proper measurement drift at the beginning of the simulation. As reference sensors, an Honeywell Miniature Inertial Measurement Unit [23] and a SSTL RIGEL-L Star Tracker [28] are taken. Their performance properties are summarized in Table 4.2. Table 4.2: Sensors performance properties.

IMU Property S M B n

Value

Scale factor Misalignment Error Bias Error ARW noise density

Property

1 170 0.005 0.005

Star Tracker Value

Relative Accuracy (1σ) Bias Error

3 5

UoM ppm µrad deg/h √ deg/ h UoM arcsec arcsec

Position and Speed Vision-based navigation systems are currently collocated between TRL 4 and TRL 6 [11], and they are still under development. A relative navigation method is required for retargeting, because guidance algorithms and hazard detection systems need informations about position and speed relatively to the Target Landing Site. Ref. [27] attests a precision better than 10 m in position and 0.2 m s−1 in velocity; in Ref. [27] a worst-case error of 2 m and 0.4 m s−1 is obtained for a near vertical landing. This kind of systems makes use of a radar or laser altimeter to estimate the altitude with which the images taken by cameras are resized to the proper scale. Since altimeters absolute error increases with the altitude, the error introduced in the model is a random error with zero-mean and variance 51

Chapter 4 proportional to the altitude. Table 4.3 shows the assumed covariance of the navigation error at 2000 m altitude, corresponding to an error of ±25 m on position and ±0.4 m s−1 on speed, with a conﬁdence level of 1σ. Table 4.3: Navigation errors covariance at 2000 m.

rx rx ry rz vx vy vz

ry

rz

vx

vy

vz

625 0 0 0 0 0 0 625 0 0 0 0 0 0 625 0 0 0 0 0 0 0.16 0 0 0 0 0 0 0.16 0 0 0 0 0 0 0.16

The error is assumed linearly decreasing with altitude up to 0 m. In the reality, lower altitudes are out of the operative range of vision-based navigation systems that are disabled during Terminal Descent [27]. The Terminal Descent phase is not covered by the simulation, so this diﬀerence can be considered admissible. Table 4.4 summarizes the assumed error covariance at ground, equivalent to an error of ±0.1 m and ±0.1 m s−1 (1σ). Table 4.4: Navigation errors covariance at ground (0 m).

rx

ry

rz

vx

vy

vz

rx 0.01 0 0 0 0 0 ry 0 0.01 0 0 0 0 rz 0 0 0.01 0 0 0 vx 0 0 0 0.01 0 0 vy 0 0 0 0 0.01 0 vz 0 0 0 0 0 0.01

4.1.3

Guidance

The guidance block receives estimations of position, speed, attitude, angular velocity and mass, and use these data, together with the coordinates of the Target Landing Site, to compute the current target attitude, angular speed and thrust magnitude, needed to reach the TLS with the desired ﬁnal states.

52

From HDA Manager rTLS rM vM qM ωM

•

GuidALG

θ ψ Tmag tToF

Guidance Control Interface

qT ωT TT

•

qT ωT qM ωM

mM

To Attitude Control

From Navigation

Retargeting Simulation

TT To Main Engine Throttle

Figure 4.4: Guidance logic diagram.

Figure 4.4 shows the structure of the guidance block. Measured states rM , vM , qM , ω M and mM are inputs for GuidALG that recompute the landing trajectory through the method expounded in chapter 3. This operation is performed: • whenever the TLS coordinates are changed by the Hazard Detection and Avoidance Manager; • anyway every 5 s, in order to cope with measure dispersions. When the spacecraft reaches an altitude of 100 m, GuidALG is disabled, and the Approach phase is completed with the last available guidance proﬁle. The GuidALG output is made by the pseudospectral proﬁles of Euler angles θ and ψ, and thrust magnitude Tmag , and by the time-of-ﬂight tToF . These values are processed by the Guidance-Control Interface block, that translates them into instantaneous target quaternions qT (t), target angular speed ω T (t) and target thrust magnitude TT (t) at every control update. First, the target attitude, in terms of target Euler angles θT (t) and ψT (t), and target thrust magnitude TT (t) are evaluated through their polynomial approximation by Eqn. (2.14). The target thrust does not require any additional calculation, and it is directly submitted to the main thrust throttle. The Euler angles time derivatives θ˙T (t) and ψ˙ T (t) are obtained through the use of the Chebyshev diﬀerentiation matrix (see section 2.2.1 at page 11) and again through the Eqn. (2.14). The roll angle φT (t) is considered constant at zero. Then, from the Eqn. (3.2) the director cosine matrix ADC is calculated from Euler angles, and it is exploited to obtain the target quaternion vector 53

Chapter 4 qT (t). From Eqn. (3.3): 1 q1 = (A23 − A32 ), 4q4 q2 = 1 (A31 − A13 ), 4q4 1 q (A12 − A21 ), 3 = 4q 4 0.5 q4 = ± 1 + A11 + A22 + A33 .

(4.11)

The target angular velocity ω T (t) is computed by solving the linear system: ˙ φ 1 − cos φ tan ψ sin φ tan ψ ωx 0 cos φ/ cos ψ − sin φ/ cos φ ωy = θ˙ , 0 sin φ cos φ ωz ψ˙

(4.12)

where the subscript T is not reported for a clearer notation. Being φ = φ˙ = 0 the system is simpliﬁed in: 0 1 − tan ψ 0 ωx 0 1/ cos ψ 0 ωy = θ˙ . 0 0 1 ωz ψ˙

4.1.4

(4.13)

Control and Actuation

Figure 4.5 shows the scheme of the Control block. Error quaternions and error angular speed are extrapolated from measures and targets. Then, a PID controller computes the theoretical control torques. These torques require to be modulated, in order to have discrete on/oﬀ commands for ACS thrusters, that they have no throttling capability. qT ωT qM ωM

Error Computation

qe ωe

PID

˜C M

PWPF Modulation

To ACS Thrusters

From Navigation and Guidance

Figure 4.5: Control logic diagram.

54

MC

Retargeting Simulation Error computation Error angular speed vector ω e is simply obtained by subtracting the target from the measure: ωe = ωM − ωT . (4.14) The rotation error is deﬁned as the rotation needed to pass from the target attitude to the measured one. For quaternions, this corresponds to:

(4.15)

qM = Rqe qT , where Rqe is: R qe

qe4 qe3 −qe2 −qe3 qe4 qe1 = qe2 −qe1 qe4 −qe1 −qe2 −qe3

qe1 qe2 . qe3 qe4

(4.16)

The system (4.15) can be rewritten isolating the error quaternion elements: qT 4 −qT 3 qT 2 qT 1 qT 3 qT 4 −qT 1 qT 2 q . (4.17) qM = −qT 2 qT 1 qT 4 qT 3 e −qT 1 −qT 2 −qT 3 qT 4

Now, the error quaternion is obtained by inverting the equation. In this case, the inverse matrix coincides with the transpose: qT 4 qT 3 −qT 2 −qT 1 −qT 3 qT 4 qT 1 −qT 2 (4.18) qe = qT 2 −qT 1 qT 4 −qT 3 qM . qT 1 qT 2 qT 3 qT 4

PID Controller ˜ C are computed from the errors, with a The theoretical control torques M generalized PID controller, using angular speed in place of derivative terms: R ˜ Cx = Kpx 2q1e q4e + Kix t 2q1e q4e dt − Kdx ωxe , M 0 Rt ˜ Cy = Kpy 2q2e q4e + Kiy 2q2e q4e dt − Kdy ωye , (4.19) M 0 Rt ˜ MCz = Kpz 2q3e q4e + Kiz 0 2q3e q4e dt − Kdz ωze ,

where the terms Kp , Kd and Ki are respectively the proportional, the derivative and the integral gains. The integral part is introduced in order to cope 55

Chapter 4 with the constant component of the thrust misalignment, and it is computed through a discrete time integrator. Error Calculation and Control blocks follow the general update rate of 20 Hz. The behaviour of a PID controller over each axis corresponds to a second order system [18] with pulsation ω and damping ξ: r Kp Kd ω= , ξ= , (4.20) I 2ωI where I is the corresponding moment of inertia over the axis. These values have been exploited to make a trade-oﬀ over gains values. The control requirements over the roll axis are less critical, and lower pulsation is allowed. Since the mass changes due to the fuel consumption, dynamics properties have been evaluated at mass mPDI = 1500 kg (the mass value at Powered Descent Initiation) and mass mdry = 790 kg. Assumed gains values and corresponding dynamics properties, are summarized in Table 4.5. Table 4.5: PID gains and equivalent dynamics properties.

Kp Roll Pitch/Yaw

Kd

Ki

1000 1000 100 1500 1500 200

m =1500 kg ω [rad/s] ξ 0.911 1.184

0.456 0.592

m =790 kg ω [rad/s] ξ 1.067 1.445

0.534 0.723

Thrusters Modulation The ACS chemical thrusters can provide only constant thrust, so the theo˜ C need to be properly modulated. A Pulse-Widthretical control torques M Pulse-Frequency (PWPF) Modulator, whose structure is showed in Figure 4.6, has been chosen. For every axis, the requested control torque by the PID controller is ﬁrstly ﬁltered by a ﬁrst order lag ﬁlter, and then passed to a Schmitt-trigger. A feedback loop subtracts the Schmitt-trigger output (that it is the modulated control torque) from the signal input. PWPF modulator oﬀers a wide range of parameters to tune and has good properties of noise rejection. Its ﬂexibility allows us to meet diﬀerent requirements through diﬀerent phases of operation, but at the same time, together with its non-linear properties, requires the tailoring of every application for each speciﬁc system, based on controller choices and system dynamics [21]. In this work, a preliminary trade-oﬀ has been made. Due to precision requirement for a pinpoint landing, together with the short duration of the 56

Retargeting Simulation

From ˜ MC PID

+

•

E −

UC Km 1 + τm s

U

•

MC

To ACS Thrusters

Figure 4.6: PWPF Modulation.

retargeting phase, minimizing phase lag has been privileged as choice criterion. Available tuning parameters are the ﬁlter gain Km , the ﬁlter time constant τm , the Schmitt-trigger on-value uon and the correspondent oﬀ-value uoff . Table 4.6 shows the actual used parameters. An output control toque umax of 40 N m is assumed available by the ACS thrusters. Table 4.6: PWPF Modulator parameters.

Parameter

Value

UoM

Km τm uon uoff umax

1 0.2 9 5 40

adim. s Nm Nm Nm

The modulation is supposed to be executed analogically by the on-board electronics, so the modulator has been built as a continuous-time Simulink block. Actuators The actuators, for both the ACS and the Propulsion System, are assumed to have a transient state much faster than the lander dynamics. Then, they are considered as ideal actuators. A minimum impulse of 20 ms for ACS thrusters is assumed.

57

Chapter 4

4.2

Results

In order to validate the model, several tests have been performed. First, a typical retargeting case is presented. Monte Carlo analysis has been exploited to verify the precision performance depending from selected landing coordinates and measures dispersions. Finally, a preliminary study on cameras Field of View has been carried out.

4.2.1

Single Landing Case

By way of example, a single typical retargeting toward the TLS at {0, 750, −1200} m from the Nominal Landing site is shown in the following pages. Figure 4.7 and 4.8 show respectively the position and the speed components during the maneuver.

X (Altitude) [m]

2000 1500 1000 500 0 0

10

20

30

40

50

60

70

80

40

50

60

70

80

40 50 Time [s]

60

70

80

Y (Downrange) [m]

1000 0 −1000 −2000 0

10

20

30

10

20

30

Z (Crossrange) [m]

0 −500 −1000 −1500 0

Figure 4.7: Single Retargeting, Position. TLS {0, 750, −1200} m.

58

Retargeting Simulation

vX (Altitude) [m/s]

0 −10 −20 −30

vZ (Crossrange) [m/s]

vY (Downrange) [m/s]

−40 0

10

20

30

40

50

60

70

80

40

50

60

70

80

40 50 Time [s]

60

70

80

40 30 20 10 0 0

10

20

30

10 0 −10 −20 −30 0

10

20

30

Figure 4.8: Single Retargeting, Speed. TLS {0, 750, −1200} m.

In Figure 4.9 the Euler angles obtained are compared to the target. The target attitude is faithfully followed. An initial phase lag is visible, subsequently corrected by following retargeting actions. This response is due to the initial discontinuity of the target angular speed, that is no explicitly considered in the trajectory computation, as can be seen in Figure 4.10, where angular velocities are displayed.

59

Chapter 4

φ [deg]

1 Target Actual

0.5 0 −0.5 0

10

20

30

40

50

60

70

80

10

20

30

40

50

60

70

80

40 50 Time [s]

60

70

80

−40

θ [deg]

−60 −80 −100 −120 0

40 ψ [deg]

20 0 −20 −40 0

10

20

30

Figure 4.9: Single Retargeting, target and actual attitude. A phase lag is visible in the early phases of the maneuver especially in the Yaw angle ψ.

60

Retargeting Simulation

ωx [rad/s]

0.04 Target Actual

0.02 0 −0.02 0

10

20

30

40

50

60

70

80

10

20

30

40

50

60

70

80

10

20

30

40 50 Time [s]

60

70

80

ωy [rad/s]

0.05 0 −0.05 −0.1 −0.15 0

ωz [rad/s]

0.05 0 −0.05 −0.1 −0.15 0

Figure 4.10: Single retargeting. Angular velocities. An initial discontinuity is visible; fast oscillations are due to the thrust misalignment disturbance torque.

870

Mass [kg]

860 850 840 830 820 810 0

10

20

30

40 50 Time [s]

60

70

80

Figure 4.11: Single retargeting. Mass trend.

61

Chapter 4 Actuators Load Figure 4.12 shows the obtained sequence of ﬁring of the ACS thrusters. The minimum commanded impulse is 20.9 ms, which satisﬁes the constraint on the minimum impulse that thrusters can supply.

YAW− YAW+ PITCH− PITCH+ ROLL− ROLL+ 0

10

20

30

40 50 Time [s]

60

70

80

Figure 4.12: Single retargeting. ACS Thruster firings sequence.

Figure 4.13 shows the trend of the commanded thrust magnitude. Small discontinuities are due to random navigation errors, that induce GuidALG to slightly change the target trajectory. Large discontinuities are a further proof of the occasional convergence of the algorithm to a suboptimal solution, subsequently ﬁxed by following retargeting checks. 2400 2200

Thrust [N]

2000 1800 1600 1400 1200 1000 0

10

20

30

40 50 Time [s]

60

70

80

Figure 4.13: Single retargeting. Main thrust magnitude trend.

62

Retargeting Simulation

4.2.2

Sensitivity to TLS

Monte Carlo simulations have been run to verify the sensitivity of the system with respect to the Target Landing Site coordinates. Every simulation considers 100 samples aﬀected by nominal navigation errors (see section 4.1.2); every simulation points toward a diﬀerent TLS. Four representative cases of diversion are presented: (a) TLS {0, 0, 0}T m (retargeting onto the Nominal Landing Site); (b) TLS {0, −1000, 0}T m (Downrange brake); (c) TLS {0, 0, −2000}T m (Crossrange diversion); (d) TLS {0, 2000, 2000}T m (composed diversion at maximum distance tested).

TLS Mean Shots 1σ 2σ 3σ

10 5

15 Z (Crossrange) [m]

Z (Crossrange) [m]

15

0 −5

10 5 0 −5 −10

−10 −15 −30

−15

−20

−10 0 10 Y (Downrange) [m]

20

30

−1030

(a) TLS {0, 0, 0}T (NLS).

−1985 −1990

TLS Mean Shots 1σ 2σ 3σ

2020

−1995 −2000

−1010 −1000 −990 Y (Downrange) [m]

−980

−970

2015 2010

2020

2030

TLS Mean Shots 1σ 2σ 3σ

2005 2000 1995

−2005

1990

−2010 −30

−1020

(b) TLS {0, −1000, 0}T .

Z (Crossrange) [m]

−1980 Z (Crossrange) [m]

TLS Mean Shots 1σ 2σ 3σ

−20

−10 0 10 Y (Downrange) [m]

(c) TLS {0, 0, −2000}T .

20

30

1970

1980

1990 2000 2010 Y (Downrange) [m]

(d) TLS {0, 2000, 2000}T .

Figure 4.14: Landing Position Sensitivity to the Target Landing Site coordinates. Initial dispersion ±25 m, ±0.4 m s−1 (1σ).

The obtained results are shown in Figure 4.14 for the position and in Figure 4.15 for the speed. The landing precision appears to be independent by the magnitude of the requested diversion. 63

Chapter 4

TLS Mean Shots 1σ 2σ 3σ

vZ (Crossrange) [m/s]

0.5 0.4 0.3 0.2 0.1 0

TLS Mean Shots 1σ 2σ 3σ

0.5 vZ (Crossrange) [m/s]

0.6

−0.1

0.4 0.3 0.2 0.1 0 −0.1 −0.2

−0.2

−0.3 −0.4

−0.2

0 0.2 0.4 vY (Downrange) [m/s]

0.6

0.8

−0.4

(a) TLS {0, 0, 0}T (NLS).

−0.2

0 0.2 0.4 vY (Downrange) [m/s]

0.6

0.8

(b) TLS {0, −1000, 0}T .

0.9

vZ (Crossrange) [m/s]

0.7 0.6 0.5 0.4 0.3 0.2 0.1

TLS Mean Shots 1σ 2σ 3σ

0.7 vZ (Crossrange) [m/s]

TLS Mean Shots 1σ 2σ 3σ

0.8

0.6 0.5 0.4 0.3 0.2 0.1 0

−0.4

−0.2

0 0.2 0.4 vY (Downrange) [m/s]

(c) TLS {0, 0, −2000}T .

0.6

0.8

−0.1

−0.4

−0.2

0 0.2 0.4 vY (Downrange) [m/s]

0.6

0.8

(d) TLS {0, 2000, 2000}T .

Figure 4.15: Landing Velocity Sensitivity to the Target Landing Site coordinates. Initial dispersion ±25 m, ±0.4 m s−1 (1σ).

64

Retargeting Simulation

4.2.3

Sensitivity to Initial Dispersion

Monte Carlo simulations have been exploited to verify the sensitivity of the system to the initial navigation dispersion. Each simulation considers 100 samples. The same targeting scenario is considered for all the simulations, with the TLS collocated at {0, 1000, 1000}T m from the Nominal Landing Site. Each simulation has a diﬀerent value of initial variance in the determination of position. Standard deviation (STD) values of 45, 25 and 10 m has been considered. Also a simulation without navigation errors has been performed, in order to determinate the impact of the navigation dispersion relatively to the accuracy obtained only from the control. TLS Mean Shots 1σ 2σ 3σ

1020 1010

1030 Z (Crossrange) [m]

Z (Crossrange) [m]

1030

1000 990 980

1020 1010

TLS Mean Shots 1σ 2σ 3σ

1000 990

960

970

980

980

990 1000 1010 1020 1030 1040 Y (Downrange) [m]

960

(a) STD 45 m (1σ).

1020 1010

TLS Mean Shots 1σ 2σ 3σ

1030

1000 990 980

980

990 1000 1010 1020 1030 1040 Y (Downrange) [m]

(b) STD 25 m (1σ).

Z (Crossrange) [m]

Z (Crossrange) [m]

1030

970

1020 1010

TLS Mean Shots 1σ 2σ 3σ

1000 990

960

970

980

990 1000 1010 1020 1030 1040 Y (Downrange) [m]

(c) STD 10 m (1σ).

980

960

970

980

990 1000 1010 1020 1030 1040 Y (Downrange) [m]

(d) No measures dispersion.

Figure 4.16: Landing Position Sensitivity to Navigation dispersions. {0, 1000, 1000}T , velocity dispersion ±0.4 m s−1 (1σ).

TLS

Figure 4.16 shows the performances obtained for the ﬁnal position; Figure 4.17 displays the results on the ﬁnal velocity. All errors due to numerical approximation in guidance are negligible (see section 3.4.1), and it can be seen that dispersion due to control is at least one order of magnitude lower than the one due to navigation errors. This proves that landing precision is mainly aﬀected by guidance errors. 65

Chapter 4

1

1 TLS Mean Shots 1σ 2σ 3σ

0.6 0.4

TLS Mean Shots 1σ 2σ 3σ

0.8 vZ (Crossrange) [m/s]

vZ (Crossrange) [m/s]

0.8

0.2 0 −0.2

0.6 0.4 0.2 0 −0.2

−0.5

0 0.5 vY (Downrange) [m/s]

1

−1

−0.5

(a) STD 45 m (1σ).

0.5

1 TLS Mean Shots 1σ 2σ 3σ

0.6 0.4

TLS Mean Shots 1σ 2σ 3σ

0.8 vZ (Crossrange) [m/s]

0.8

0.2 0 −0.2

0.6 0.4 0.2 0 −0.2

−1

−0.5

0 vY (Downrange) [m/s]

(c) STD 10 m (1σ).

0.5

1

−1

−0.5

0 vY (Downrange) [m/s]

0.5

(d) No measures dispersion.

Figure 4.17: Landing Velocity Sensitivity to Navigation dispersions. {0, 1000, 1000}T , velocity dispersion ±0.4 m s−1 (1σ).

66

1

(b) STD 25 m (1σ).

1

vZ (Crossrange) [m/s]

0 vY (Downrange) [m/s]

TLS

Retargeting Simulation

4.2.4

Cameras Field of View

A preliminary study on cameras Field of View (FoV) for navigation and/or hazard detection has been performed. The smooth attitude proﬁle imposed by the trajectory polynomial approximation guarantees the absence of sudden maneuvers, allowing a continuous landmarks tracking. A single camera, pointing towards the roll axis (xB ) has been assumed. By way of example, Figure 4.18 shows how the intersection between the camera line of sight and the ground varies during a representative retargeting (TLS at {0, 750, −1200}T m). The camera pointing depicted into the graph is taken every 0.5 s, an update rate lower than the actual camera frame rate, that is considered between 10 and 20 Hz [11]. Tracked landmarks pass into the FoV in a continuous way, allowing relative navigation. Camera Target Ground Track TLS NLS

800 600

z (Crossrange) [m]

400 200 0

START

−200 −400 −600 −800 −1000 −1200 −1500

−1000

−500 0 Y (Downrange) [m]

500

1000

Figure 4.18: Camera pointing during retargeting. Blue crosses represent the intersection between the camera line of sight and the ground, taken every 0.5 s during the landing.

The nominal landing trajectory is designed in order to have the necessary time for hazard mapping before the ﬁrst retargeting. But if a second retargeting is required, during the ﬁrst maneuver the TLS is required to be in sight for the time necessary to update the hazard map. Figure 4.19a shows the trend of the angle between the line of sight and the TLS direction during the same maneuver of Figure 4.18. Since the FoV angle considered is between 50 and 70◦ [11], the LS visibility can be considered lost when this TLS-Sightline angle is greater than 25-35◦ , depending by the actual HDA system architecture. It can be seen that the oscillatory movement implied by the diversion causes this loss just after the diversion start. However, the visibility is recovered at lower altitudes, enabling the recreation of the hazard map. This 67

Chapter 4 behaviour is obtained in all the analysed cases. For minor retargetings, the oscillation magnitude is smaller, and the LS visibility can be maintained along the the maneuver, as showed for a diversion of {0, 150, 150}T m in Figure 4.19b. Diversions up to ±2000 m on both the horizontal axes have been tested. In diversions above 1500 m, greater oscillations are counterbalanced by greater time-of-ﬂight, that guarantees more time for a second retargeting, as displayed in Figure 4.19c. 30

100

Angle [deg]

Angle [deg]

80 60 40 30 deg 20 0 0

30 deg

25 20 15 10 5

10

20

30

40 Time [s]

50

60

70

80

0 0

(a) TLS at {0, 750, −1200}T m.

10

20

30 40 Time [s]

50

60

70

(b) TLS at {0, 150, 150}T m.

120

Angle [deg]

100 80 60 40

30 deg

20 0 0

10

20

30

40 50 Time [s]

60

70

80

90

(c) TLS at {0, 1700, 1600}T m.

Figure 4.19: Angle between the camera line of sight and the TLS direction. At great angles the landing site might be out from camera FoV.

If a double retargeting is expected, a change in the retargeting strategy, in order to guarantee the necessary visibility, could be considered. One possibility is a two-step retargeting similar to that suggested in [17]: the ﬁrst diversion is commanded at high altitude (2000 m or above if it is possible), pointing towards the vertical onto selected TLS at an intermediate altitude (400-600 m). In this way the system is able to perform, at the end of this maneuver, a short vertical descent at constant speed, during which the hazard map can be updated with more resolution. Then, if requested, a second diversion until the TD phase and the touchdown can be commanded. 68

Retargeting Simulation If not needed, the descent can be completed with a fuel optimal vertical descent, whose solution is known in close form [22]. Another possible choice is a modiﬁcation in the hazard detection system architecture, such as the use of multiple cameras, or the use of gimbals, in order to increase the Field of View [27].

69

Chapter

5

Conclusions The purpose of this work was the development of a retargeting algorithm, capable of updating and correcting a landing trajectory almost to the touchdown in a planetary powered descent. A classical polynomial approach has been extended, in order to improve ﬂexibility in the landing site choice, and to consider additional non linear constraints during the descent. When a retargeting is ordered, a polynomial approximation of the optimal solution is found by a compass search algorithm. Divert capabilities, computational and precision performances have been tested through Monte Carlo methods. The functionality and the robustness of the algorithm have been tested by applying it in a simulation of a complete retargeting. The guidance scheme has been coupled with an attitude controller, and perturbed states have been exploited in order to emulate navigation system errors. In order to identify possible sources of errors in placing the spacecraft on target, a dispersion analysis has been performed. It has been observed that navigation errors play a major role in determining the accuracy at touchdown.

5.1

Considerations

A loss of divert capability in braking, located at mid-low altitude, has been observed (see section 3.4.2). In fact, under an altitude of 1000 m, the attainable landing area is no longer centered on NLS, but moves onward along the Downrange direction, reducing the capabilities of a shorter landing. This phenomenon arises in phases of the nominal trajectory in which the main engine is requested to supply an almost maximum thrust. This implies a decreased manoeuvrability in all the contexts in which additional thrust is required (as in braking). In order to regain divert capability, some available 71

Chapter 5 thrust could be deserved for retargeting, or a diﬀerent retargeting strategy could be considered, as the double maneuver scheme discussed in section 4.2.4. Some additional considerations about the Terminal Descent phase have to be made. Since the proposed guidance algorithm does not take explicitly into account rotational velocities, the earlier instants of the TD phase must be dedicated to cancel the residual angular speed in the shortest possible time. In fact, during Terminal Descent the spacecraft is required to have a constant vertical attitude and a constant vertical speed until the touchdown.

5.2

Future Developments

The retargeting problem presents several topics deserving of further investigations: 1. Retargeting Strategy In order to satisfy visibility constraints, or to improve retargeting capability, a diﬀerent retargeting strategy could be considered, as the previously mentioned two step maneuver. In that case, a trade-oﬀ on intermediate boundaries conditions would be necessary, in order to guarantee good performance with respect to both fuel consumption and ﬂexibility of landing site selection. 2. Trajectory Formulation Terms of higher degree can be added to polynomial formulation, in order to satisfy additional constraints or improve the algorithm ﬂexibility, giving more complexity to the attainable trajectory. In contrast, new optimization variables could increase the algorithm computational cost. It is known that the most general optimal thrust proﬁle is discontinuous with a max-min-max structure. The introduction of multiple polynomial domains, linked by continuity conditions, could be considered to improve the optimality of the solution. 3. Optimization Algorithm This work has shown that the proposed approach allows us to reach a wide divert capability with a simple optimization algorithm. Anyway, this method can occasionally converge towards suboptimal solutions and it is susceptible of improvement: 72

Conclusions • An intense trade-oﬀ activity should be conducted to determine the best values of the algorithm parameters: maximum iteration values, initial mesh size, minimum allowed mesh size, growth/decrement ratios of mesh size, type of generating set. All these parameters aﬀect computation time and convergence properties. • A proper ﬁrst-guess solution reduces the number of iterations required for convergence. Both achieved optimality and computational speed are aﬀected by the eﬀectiveness of the initial guess. An option could be the choice of a classical Apollo-like polynomial trajectory. • Others direct search methods could be investigated, in order to improve optimality, convergence speed and feasibility region. 4. Code Efficiency R code, All the algorithm and simulations have been written in Matlab which is a fourth-generation programming language, and tested onto a desktop PC. A more realistic estimate of the computational performances achieved could be obtained by testing algorithms on real hardware. Finally, further developments cannot be carried out without considering the actual performances and requirements of Navigation and Hazard Detection systems. The natural prosecution of this work is the integration of these two subsystem with Guidance and Control module, from numerical simulations up to hardware-in-the-loop testing.

73

Bibliography

[1] B. Acikmeşe and S. R. Ploen. “A Powered Descent Guidance Algorithm for Mars Pinpoint Landing”. In: AIAA Guidance, Navigation, and Control Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2005. [2] B. Acikmeşe and S. R. Ploen. “Convex Programming Approach to Powered Descent Guidance for Mars Landing”. In: Journal of Guidance, Control and Dynamics 30.5 (2007). [3] M. Alger et al. “Pinpoint Lunar Landing Navigation using Crater Detection and Matching: Design and Laboratory Validation”. In: AIAA Guidance, Navigation, and Control Conference. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2012. [4] F. V. Bennett. “Lunar descent and ascent trajectories”. In: 8th Aerospace Sciences Meeting. Aerospace Sciences Meetings. American Institute of Aeronautics and Astronautics. New York, 1970. [5] John T. Betts. “Survey of Numerical Methods for Trajectory Optimization”. In: Journal of Guidance, Control, and Dynamics 21.2 (1998), pp. 193–207. [6] C. Canuto et al. Spectral methods in fluid dynamics. New York: Springer, 1988. [7] J. D. Carpenter et al. “Scientiﬁc preparations for lunar exploration with the European Lunar Lander”. In: Planetary and Space Science 74.1 (2012), pp. 208–223. [8] ESA Ministerial Council 2012 on ESA website. 2013. url: http:// www.esa.int/About_Us/Ministerial_Council_2012. 75

[9] ESA Lunar Lander on ESA website. 2013. url: http://www.esa.int/ Our_Activities/Human_Spaceflight/Lunar_Lander. [10] F. Fahroo and I. M. Ross. “Direct Trajectory Optimization by a Chebyshev Pseudospectral Method”. In: Journal of Guidance, Control and Dynamics 25.1 (2002). [11] G. Flandin et al. “Maturing Vision Based Navigation Solutions to Space Exploration”. In: AIAA Guidance, Navigation, and Control Conference. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2010. [12] B. Gardini. “OUTLOOK after the CM 2012”. In: EC Workshop on Space Science and Exploration. European Commission. Madrid, 2013. url: http://ec.europa.eu/enterprise/newsroom/cf/itemdetail. cfm?item_id=6304&lang=en. [13] A. M. Hawkins. “Constrained Trajectory Optimization of a Soft Lunar Landing From a Parking orbit”. M.Sc. 2005. [14] N. Henn et al. “Chemical Propulsion Challenges for ESA’s Exploration Missions and Automatic Return Vehicle”. In: 47th AIAA/ASME/SAE/ ASEE Joint Propulsion Conference & Exhibit. Joint Propulsion Conferences. American Institute of Aeronautics and Astronautics. 2011. [15] S. Hobbs, D. Rosa, and J. Delaune. “Guidance and Control system design for Lunar Descent and Landing”. In: AIAA Guidance, Navigation, and Control Conference. Guidance, Navigation, and Control and Co-located Conferences. 2010. [16] A. E. Johnson et al. “Analysis of On-Board Hazard Detection and Avoidance for Safe Lunar Landing”. In: Aerospace Conference, 2008 IEEE. ID: 1. Institute of Electrical and Electronics Engineers. 2008, pp. 1–9. [17] M. Johnson. “A Parameterized Approach to the Design of Lunar Lander Attitude Controllers”. In: AIAA Guidance, Navigation, and Control Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2006. [18] M. H. Kaplan. Modern Spacecraft Dynamics & Control. New York: John Wiley & Sons, 1976. [19] A. R. Klumpp. “Apollo lunar descent guidance”. In: Automatica 10.2 (1974), pp. 133–146.

76

[20] T. Kolda, R. Lewis, and V. Torczon. “Optimization by Direct Search: New Perspectives on Some Classical and Modern Methods”. In: SIAM Review 45.3 (2003), pp. 385–482. [21] T. D. Krøvel. “Optimal Tuning of PWPF Modulator for Attitude Control”. M.Sc. 2005. [22] J.S. Meditch. “On the problem of optimal thrust programming for a lunar soft landing”. In: Automatic Control, IEEE Transactions on 9.4 (1964), pp. 477–484. [23] Miniature Inertial Measurement Unit. Datasheet. Honeywell. 2009. url: http://www51.honeywell.com/aero/common/documents/ myaerospacecatalog-documents/Space-documents/Miniature_ Inertial_MeasurementUnit.pdf. [24] B. Parreira et al. “Consolidated Performance Assessment of Hazard Avoidance Techniques for Vision Based Landing”. In: AIAA Guidance, Navigation and Control Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. Hilton Head, South Carolina, 2007. [25] R. Peyret. Spectral Methods for Incompressible Viscous Flow. Applied Mathematical Sciences. New York: Springer, 2002. [26] A. Pradier et al. “The ESA Lunar Lander Mission”. In: AIAA SPACE 2011 Conference & Exposition. SPACE Conferences & Exposition. American Institute of Aeronautics and Astronautics. 2011. [27] J. Riedel et al. “Optical Navigation Plan and Strategy for the Lunar Lander Altair; OpNav for Lunar and other Crewed and Robotic Exploration Applications”. In: AIAA Guidance, Navigation, and Control Conference. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2010. [28] RIGEL-L Star Tracker. Datasheet. Surrey Satellite Technology Ltd. 2013. url: http://www.sstl.co.uk/Products/Subsystems/ Flying-Soon/Rigel-L-Star-Tracker. [29] L. Singh and S. Lim. “On Lunar On-Orbit Vision-Based Navigation: Terrain Mapping, Feature Tracking Driven EKF”. In: AIAA Guidance, Navigation and Control Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2008. [30] Soyuz from Guyana Space Center User’s Manual. Arianespace. 2012. url: http://www.arianespace.com/launch-services-soyuz/ Soyuz-Users-Manual-March-2012.pdf. 77

[31] D. H. Titterton and J. L. Weston. Strapdown inertial navigation technology. 2nd. Progress in astronautics and aeronautics. Reston, VA; Steveneage, UK: American Institute of Aeronautic, Astronautic; Institution of Electrical, and Electronis Engineers, 2004. [32] U. Topcu, J. Casoliva, and K. Mease. “Fuel Eﬃcient Powered Descent Guidance for Mars Landing”. In: AIAA Guidance, Navigation, and Control Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2005. [33] L. N. Trefethen. Spectral methods in Matlab. Philadelphia, 2000. [34] D. Unsal and K. Demirbas. “Estimation of deterministic and stochastic IMU error parameters”. In: Position Location and Navigation Symposium (PLANS), 2012 IEEE/ION. Institute of Electrical and Electronics Engineers. 2012. [35] E. Wong, J. Masciarelli, and G. Singh. “Autonomous Guidance and Control Design for Hazard Avoidance and Safe Landing on Mars”. In: AIAA Atmospheric Flight Mechanics Conference and Exhibit. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics. 2002.

78

Copyright © 2019 PROPERTIBAZAR.COM. All rights reserved.