Study and Industrialization of Optimization Methods for Low-Thrust Orbital Maneuvers

(1)

IN

DEGREE PROJECT MATHEMATICS, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2016,

Study and Industrialization of Optimization Methods for Low- Thrust Orbital Maneuvers

ROMAIN BOYER

KTH ROYAL INSTITUTE OF TECHNOLOGY SCHOOL OF ENGINEERING SCIENCES

(2)

(3)

Study and Industrialization of Optimization Methods for Low-Thrust Orbital Maneuvers

R O M A I N B O Y E R

Master’s Thesis in Systems Engineering (30 ECTS credits) Degree Programme in Vehicle Engineering (120 credits) Royal Institute of Technology year 2016 Supervisor at Thales Group: Jérémie Labroquère Supervisor at KTH: Per Engvist Supervisor at ISAE: Stéphanie Lizy-Destrez Examiner at KTH: Per Engvist

TRITA-MAT-E 2016:72 ISRN-KTH/MAT/E--16/72--SE

Royal Institute of Technology School of Engineering Sciences KTH SCI SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(4)

(5)

Acknowledgements

I would like to use this opportunity to thank all people for their welcome inside the friendly but working atmosphere of the Flight Dynamics department.

I want to address special thanks to Jérémie Labroquère who gave his time to supervise me during these six months. His advice, explanations and comments have been very helpful.

I would like also to thank all people who answered my questions, who helped me in any way to master a bug with the dierent software I had the opportunity to use during the internship, or who gave me advice about orbital maneuver and optimal control: Thomas Ro- drigues, Emmanuel Bignon, Pierre Mercier and Johann Tournebize.

Finally, I am grateful to the two professors Stéphanie Lizy-Destrez from ISAE and Per Enqvist from KTH for accepting to be supervisors and examiners of this degree project.

(6)

(7)

CONTENTS CONTENTS

6 Encountered problems 63 7 Conclusion 64 Appendices 65 A Code in C++ using solver Nlopt (SLSQP algorithm) 65 B Theoretical solution of the Brachistochrone problem 66 C Output of Ipopt solver with the Brachistochrone problem (single shooting) 67 D Test problems for orbital transfers 69 D.1 Brachistochrone problem . . . 69

D.2 Van der Pol oscillator . . . 69

D.3 Maximum-Energy Orbit Transfer . . . 69

D.4 Maximum-Radius Orbit Transfer with Solar Sail . . . 69

D.5 Planet to planet transfer: Earth-Mars . . . 70

D.6 Earth-to-Mars problem with Nuclear-electric propulsion systems . . . 70

D.7 Earth-to-Mars problem with Solar-electric propulsion systems . . . 71

D.8 Earth-to-Mars problem with Solar sail spacecraft . . . 71 E Results for a high-thrust transfer from the paper "Low-Thust, High-Accuracy

Trajectory Optimization" 72

F Results for a medium-thrust transfer from the paper "Low-Thust, High-

Accuracy Trajectory Optimization" 73

G Output of Ipopt solver with a low-thrust problem (change of SMA and

single shooting method) 74

References 75

(9)

LIST OF FIGURES LIST OF FIGURES

List of Figures

1 Thales Services into Thales Group . . . 6

2 Thales Services - Organization . . . 6

3 Organization of LOM . . . 7

4 Planning . . . 9

5 Keplerian elements . . . 11

6 Sketch of an optimal control problem . . . 13

7 Sketch of the single shooting method . . . 17

8 Sketch of the multiple shooting method . . . 18

9 Brachistochrone solution with Nlopt (Collocation) . . . 29

10 Brachistochrone solution with Ipopt (Single shooting) . . . 30

11 Brachistochrone solution with Ipopt (Multiple shooting) . . . 31

12 General view of the software architecture: Package diagram. . . 35

13 General view of the software architecture: Class diagram. . . 37

14 The "direct method" package. . . 38

15 Optimal radial thrust for the time-optimal problem . . . 41

16 Optimal tangential thrust for the time-optimal problem . . . 41

17 Evolution of the main parameters . . . 42

18 Transfer trajectory . . . 42

19 Optimal radial thrust of the time-optimal problem . . . 43

20 Optimal tangential thrust of the time-optimal problem . . . 43

29 Optimal vertical thrust of the time-optimal problem . . . 52

47 Theoretical solution of the Brachistochrone problem . . . 66

48 Optimal solution for a time-optimal problem . . . 72

49 Optimal solution for a time-optimal problem . . . 73

(10)

List of Tables

1 Avantages and disavantages of direct methods . . . 15

2 NLP Solvers (part 1/3) . . . 20

3 NLP Solvers (part 2/3) . . . 21

4 NLP solvers (part 3/3) . . . 22

5 NLP characteristics according to the selected direct method for the Brachis- tochrone problem . . . 32

6 Orbits characteristics . . . 40

8 Physical Constants . . . 50

(11)

1 INTRODUCTION

1 Introduction

1.1 Presentation of the company 1.1.1 Thales Group and Thales Services

Thales Group, present in 56 countries and employing around 67000 people, is one of the lead- ers in critical information systems in the Aerospace, Transportation and Defense and Security markets. The group is mostly known for its broad experience of the development and synergy of state-of-art technologies both for the civil and military markets, and also for its Research and Development activities which bring together a multinational group of around 25000 en- gineers and researchers.

Present in civil and military eld, Thales Group is organized into 6 divisions: Secure Commu- nications and Information Systems, Avionics, Land and Air Systems, Space, Defense Mission Systems and Ground Transportation Systems. Each division is divided into Strategic Business Units (SBU), and the latters into Strategic Business Lines (SBL). The activities of Thales Services are part of the Secure Communications and Information Systems division.

Figure 1: Thales Services into Thales Group

The following gures show more precisely the organization of Thales Services:

Figure 2: Thales Services - Organization

(12)

Figure 3: Organization of LOM 1.1.2 Flight dynamics department

Inside Thales Service, the departments are divided into a matrix organization considering the geographic area and the domain of activity.

The internship has been done in the Flight Dynamics Department in Toulouse, which is part of Division South-West, inside the "Logiciel Orienté Machine" Unit (LOM).

The Flight Dynamics department, as the name indicates, provides services in the area of space mechanics, mainly to serve spacecraft operation's needs. In particular, the department provides two types of services:

• Technical support in the client's facilities, concerning the areas of mission analysis, technical studies

• Software project development directly at Thales site, which includes in particular the development of software tools such as:

- CMS (Space Mechanic Center): Space Mechanics library used by CNES (for control operations) and developed in FORTRAN. It includes the low-level libraries and high level services as maneuvers computations, guidance, orbit determination (OD), commands generation and so forth

- PATRIUS: Analogue to CMS but developed in Java - ZOOM : Orbit Determination toolkit developed for CNES - GOTLIB/PHRLIB (C++): Attitude Guidance Libraries

- STELA: Long-term propagation tool developed for end-of-life activities

- GALILEO: Guidance library for the next generation European navigation system - ASPROSH: Commands Validation Library

- COSMEUS (JACK, JMAN, JPOD, JPAD): Software suite for space mechanics - POD/PAD for Gokturk mission: Precise Orbit Altitude Determination Library

(13)

1.2 Internship description 1 INTRODUCTION

1.2 Internship description 1.2.1 Context

This internship extends a Master Thesis, which was written by Hanane Ait-Lakbir in 2015 entitled "Study and Industrialization of Computational Methods for Orbital Maneuver" [15].

Hanane's Master Thesis was a part of the development software suite COSMEUS (Complex Space Mechanics Utilities and Software), a java library related to space mechanics which is composed of the following modules:

• JACK (Java Aerospace Capitalization Kit): low-level library for space mechanics (orbit propagation, force models, orbital event management, orbital parameter computation...)

• JUDE (Java Utilities for Development): library for exception, message, logs management

• JSCAR (Java Software for Comparison, Analysis and Reporting): library for tests

• JMAN (Java Maneuver): library for maneuver computation

• POD (Java Precise Orbit Determination) and JPAD (Java Precise Attitude Determina- tion): generic library for orbit and attitude determination

• Other modules: JERY for Java Earth Reentry: library for atmospheric reentry. JDEB (Java Debris) for debris trajectory propagation and assessment of collision risks.

Hanane's Master Thesis completed JMAN and JACK with some functionnalities in order to integrate electrically propelled maneuvers (ie low-thrust maneuvers). More precisely, Hanane used indirect optimization methods to compute optimal maneuvers. On the contrary, this Master Thesis deals with direct methods to optimize low-thrust orbital tranfers. It is quite important to know how to use thrust during a space mission in order to optimize one parameter. For example, for a space company, one could want to minimize the fuel consumption, with respect to some mission constraints.

1.2.2 Flight trajectories computation issue

In this part, the ight trajectories computation issue is explained. First, orbital maneuvers can be divided into four dierent types:

• Station-keeping: the aim is to maintain the satellite on its orbit because some pertubations modify the trajectory

• Orbit raising and early orbit phase: trajectory between two dierent orbits

• Debris avoidance: to avoid collisions which could damage the satellite

• Interplanatery transfers: between two planets, typically a mission to go to Mars Only orbit raising maneuvers are studied in this report. There could be trajectories between two orbits with changes in SMA (semi-major axis), changes in eccentricity, changes in inclination or transfers from GTO (Geostationary Transfer Orbit) to GEO (Geosynchonous Orbit).

This kind of maneuvers is performed with two dierent types of propulsion: chemical or electric propulsion. Nowadays, the electric propulsion (ie low-thrust propulsion) is becoming more and more important because it allows to consume less propellant. Thus, lighter satellites can be built which save some money for the aerospace industry. However, the main disavantage is the duration of the maneuvers which is signicantly higher. For example, a spiral transfer, between a parking orbit at 200 km to GEO, with a Isp equal to 310 s and a constant ratio _m^T equal to 0.01 m.s⁻², can last ve days instead of ve hours with a Hohmann transfer.

(14)

During a space mission, some parameters can be optimized such as the fuel consumption, the radius of the nal orbit or the duration of the maneuver, with respect to some space mission constraints. In the case of chemical propulsion, it is quite easy to compute the optimal maneuvers because a maneuver can be modeled as a sequence of impulsive thrusts. Indeed, the duration of thrust arcs is really short compared to the mission duration. However, for low-thrust propulsion (ie electric propulsion), optimization problems are less straightforward to solve. The control variables (thrust) become continuous functions and make the problem harder to solve. Continous optimization tools are then needed.

1.2.3 Internship purposes

The aim is to assess the possibility to use direct optimization methods in order to solve low- thrust trajectory problems. Based on some test cases, the applicability to compute optimal orbital transfers with low-thrust propulsion and the potential replacement of indirect methods will be detailed. The focused maneuvers are orbit raising maneuvers, for example with changes in eccentricy or SMA. Either the fuel consumption or the duration of the trajectory can be optimized with respect to some mission constraints.

1.2.4 Planning

The initial planning can be observed below:

Figure 4: Planning

This planning has almost not changed during the internship, except for the part "get used to Jman" which lasted only two weeks. As you can see, most of the time is dedicated to the implementation and the validation of a software.

(15)

2 MODELLING AND OPTIMIZATION TECHNIQUES

2 Modelling and optimization techniques

In this section, the two-body problem is rst described. Then, the other pertubations are detailed to obtain a better approximation of what could be the motion of a satellite in the near- Earth space. Lastly, the optimization theory is briey explained with the general formulation of an optimal control problem and dierent methods of resolution.

2.1 Satellite ight trajectory modelling 2.1.1 Two-body problem and Keplerian motion

The gravitational two-body problem concerns the motion of two point particles that interact only with each other, due to gravity. Any other inuence is neglected. In this master thesis, the studied two-body problem is given by the Earth and the spacecraft. Indeed, only the raising orbit maneuvers around the Earth are considered as it was said before. For example, interplanatery transfers could have included other bodies like the Sun or Mars. This model allows to approximate the motion of a spacecraft around the Earth by considering only one force which is gravity and using the Kepler's law.

Let consider the Earth, which is assumed to be xed, and a satellite. The following parameters can be introduced in a cartesian frame:

r = (x, y, z) position vector v = ˙r = ( ˙x, ˙y, ˙z) velocity vector a = ˙v = (¨x, ¨y, ¨z) acceleration vector

As it was said before, only one force (gravity) is taken into account in this two-body problem.

Using the Newton's law of gravity, the gravitation eld can be expressed as:

F = −G.M.m

r³ .r = −µ.m r³ .r

where r corresponds to the distance between the satellite and the Earth center, G the univer- sal constant of gravitation and µ the gravitational parameter.

From there, the equation of motion becomes:

¨ r + µ

r³.r = 0

The Kepler's laws can now be stated in the case of a satellite around the Earth. The proofs of the three laws are not written because it does not seem to be relevant here:

• The orbit of a spacecraft is an ellipse with the Earth at one of the two foci.

• A line segment joining a satellite and the Earth sweeps out equal areas during equal intervals of time.

• The square of the orbital period of a satellite is proportional to the cube of the semi- major axis of its orbit.

Instead of the regular cartesian elements, it may be very useful to use orbital elements to identify a specic orbit. Traditionnally, these elements are the six Keplerian elements:

• Eccentricity e

• Semi-major axis a

• Inclination i (angle between the reference plane and the orbital plane)

(16)

• True anomaly θ

• Right ascension of the ascending node Ω (angle between the vernal axis and the ascending node axis)

• Argument of the perigee ω (angle between the direction of the perigee and the ascending node axis).

These elements are illustrated on a specic orbit with the next gure:

Figure 5: Keplerian elements

Nevertheless, the Keplerian parameters are not the most suitable for electric transfers because they are not dened in some cases, for example if e tends to zero. On the contrary, the Equinoctial parameters are very useful to model a low-thrust problem (cf section 5). Indeed, the equinoctial parameters often help to avoid some numerical problems such as a division by zero or scaling problems. These parameters are p, which is the semi-parameter, f and g which describe the eccentricity of the orbit, h and k which describe the inclination, L which is the true longitude. They can be related to the keplerian parameters with the following equations:

p = a(1 − e²) f = e cos(ω + Ω) g = e sin(ω + Ω) h = tan(i

2)cos(Ω) k = tan(i

2)sin(Ω) L = ω + Ω + θ

(17)

2.2 Optimal Control Theory 2 MODELLING AND OPTIMIZATION TECHNIQUES

The conversion towards the cartesian elements is ensured with the following equations:

α² = h²− k² s² = 1 + h²+ k²

w = 1 + f.cos(L) + g.sin(L) r = p

w x = r

s²(cos(L) + α².cos(L) + 2.h.k.sin(L)) y = r

s²(sin(L) − α².sin(L) + 2.h.k.cos(L)) z = r

s²(h.sin(L) − k.cos(L)) 2.1.2 Perturbations

The gravitation eld is not the only force which aects the satellite. Consequenltly, the keplerian motion can only give an approximation of the trajectory. The main perturbations of the path for a satellite in the near-Earth space can be summarized as below:

• Earth shape (not perfectly spherical): the eects of the Earth atness change the expression of the gravitation eld.

• Other bodies' inuences, mainly the Sun and the Moon

• Atmosphere (drag)

• Solar radiation: apparition of a pressure force due to the exchange of quantity of motion between the photons and the spacecraft.

2.2 Optimal Control Theory

A trajectory optimization problem can generally be seen as an optimal control problem and solved using indirect or direct methods. In this master thesis, only direct methods are used.

2.2.1 General formulation

The aim of an optimal control problem is to determine a control law u(t), related to a nite state x(t) with some dierential equations corresponding to state dynamics:

˙

x(t) = f (t, x(t), u(t))

with boundary conditions:

x(t₀) = x₀ initial condition

ψ_i(xtf) ≤ 0 inequality f inal condition ψ_e(x_t_f) = 0 equality f inal condition

path constraints:

K_i(t, x(t), u(t)) ≤ 0 inequality path constraint Ke(t, x(t), u(t)) = 0 equality path constraint

(18)

and box constraints:

ulower≤ u(t) ≤ u_upper x_lower≤ x(t) ≤ x_upper

in order to optimize a cost function:

Z tf

t0

L(x(t), u(t)) dt + φ(t_f, x_f)

L is a C¹ function called the cost rate, φ a C⁰ function called the nal cost. The general formulation easily follows for a minimization:

min Z _t_f

t0

L(x(t), u(t)) dt + φ(t_f, x_f) cost f unction st x(t) = f (t, x(t), u(t))˙ dynamics

x(t0) = x0 initial condition

ψ_i(x_t_f) ≤ 0 inequality f inal condition ψe(xtf) = 0 equality f inal condition

K_i(t, x(t), u(t)) ≤ 0 inequality path constraint Ke(t, x(t), u(t)) = 0 equality path constraint ulower≤ u(t) ≤ u_upper box constraint xlower≤ x(t) ≤ x_upper box constraint

The following gure illustrates the previous general form of an optimal control problem.

Figure 6: Sketch of an optimal control problem For a maximization, it would have been "min −

Rtf

t0 L(x(t), u(t)) dt + φ(t_f, x_f)" instead of "min

Rtf

t0 L(x(t), u(t)) dt + φ(t_f, x_f)" in the cost function.

To be more concrete, in the case of trajectory optimization, the control law u(t) corresponds to the thrust whereas x(t) could include both position and velocity of the satellite. Moreover, the dierential equations are traditionnally the equations of motion and the cost function could be the mission duration. All this notions related to low-thrust transfers are described later.

(19)

2.2.2 Solving methods of an optimal control problem There are four types of methods:

• Dynamic programming

• HJBE (The HamiltonJacobiBellman Equation)

• Indirect methods (Calculus of variations and PMP)

• Direct methods (Sequential methods and simultaneous methods)

The two last methods are traditionally used for low-thrust trajectory optimization.

2.2.2.1 Dynamic Programming

This method provides the optimal closed-loop control law. With dynamic programming, you can obtain a global solution. However, it is barely used for the computation of orbital maneuvers because the computation of the maneuver involves the reaction of perturbations.

2.2.2.2 HJBE

This method corresponds to a continuous-time dynamic programming and provides satisfactory properties such as sucient conditions for optimality and optimal solutions in feedback form. However, there are also several drawbacks. Mainly, the HJB equations is very hard to solve numerically and analytic solutions can only be obtained in a few cases.

2.2.2.3 Indirect methods

They are based on PMP (Pontryagin's Maximum Principle) and Calculus of Variations. They were used in Hanane's Master Thesis to compute optimal solutions for low-thrust transfers.

These methods allow to reach high accuracies and smooth solutions. Nevertheless, some disavantages can be stated. Indirect methods are quite complex to use. It is quite dicult to nd a global solution. Lastly, the main problems remain the convergence of the shooting methods and the sensitivity of the initial guess. All these diculties explain the choice to use direct methods in this Master Thesis.

Let us describe more precisely indirect methods by showing how to solve PMP in practice.

• Dene the Hamiltonian function:

H(t, x, u, λ) = L(t, u(t)) + λ^Tf (t, x(t), u(t))

• Pointwise minimization of H with respect to u:

˜

µ(x, λ) = arg min

u H(x, u, λ) The optimal control is expressed by:

u^∗(t) = ˜µ(x(t), λ(t))

• Solve TPBVP (Two Point Boundary Value Problem):

λ = −˙ ∂H

∂x(t, x(t), ˜µ(x(t), λ(t)), λ(t)),

˙

x(t) = f (x, ˜µ(x(t), λ(t))),

(20)

with boundary conditions

λ(t_f) − ∂φ

∂x(x(t_f)) ⊥ S_f, x(t_i) = x_i, x(t_f) ∈ S_f.

TPBVP can be numerically solved with shooting techniques (single or multiple shooting).

2.2.2.4 Direct methods

A direct method tranforms an optimal control problem into a non-linear problem (NLP), which can normally be solved by an optimization solver. The following table sums up the characteristics of a direct method, which can be compared with a indirect method in the previous paragraph.

Advantages Disavantages

Convergence of the algorithm Existence of local mimimum solutions Choice of the initial guess not too crucial Time consuming if large dimensions

Quite robust algorithm Relative low precision Table 1: Avantages and disavantages of direct methods

2.2.2.5 Hybrid methods

Hanane Ait-Lakbir has used homotopy continuation methods to solve the diculties related to the indirect methods. One of the main problems was the convergence of the algorithm. A possibility could be to use hybrid methods (both direct and indirect methods). More precisely, one can apply PMP to get a boundary-value problem. Then, direct methods can be used to initialize a multiple shooting method. Direct methods allows to obtain a good initial guess for PMP methods. In other words, they provide a point not too far away from the global solution to avoid convergence problems. The advantage of this approach is to provide an initial guess close to the global solution. In this way, it is possible to get a better accuracy than direct methods and avoid the convergence problems of indirect methods.

• First step: PMP to get TPBVP

• Second step: Resolution of the initial optimal control problem thanks to a direct method.

• Third step: Initialization of the multiple shooting with the point given by the direct method.

• Fourth step: Resolution of TPBVP with multiple shooting.

2.2.3 Direct methods applied to optimal control problem

In this Master Thesis, only direct methods are performed to solve trajectory optimization problems with low-thrust propulsion. This choice has been motivated by the main avantages of this kind of methods which are described in table 1. The direct methods can be divided into two dierent categories:

• Sequential methods: only the control is discretized. The collocation is likely to be the most famous sequential method.

(21)

• Simultaneous methods: both the state vector and the control are discretized. It often leads to a direct shooting, which can be single or multiple.

The aim is to get a non-linear problem (NLP) by discretizing both the state vector and the control in the initial optimal control problem (simultaneous approach). The rst step is to discretize the time interval t ∈ [0, tf] into N − 1 equal subintervals (ie N points). The step length is given by h = _{N −1}^t^f so that tj+1− t_j = h, j = 0, . . . , N − 1. Both the state vector and the control are discretized by:

x(t) = x_j, u(t) = u_j, where jh ≤ t < (j + 1)h.

It leads to N points for the state vector and N-1 points for the control because the control uis not considered at the nal point. Let dimu be equal to the dimension of the control u, dim_x the dimension of the statevector x. Both matrices X and U, which respectively contain all the state vector coecients and all the control coecients for all the steps, are introduced:

X =





x_0,0 ... x_{0,N −1}

... ... ...

x_dim_x−1,0 ... x_dim_x−1,N −1



U =





u_0,0 ... u_{0,N −2}

... ... ...

u_dim_u−1,0 ... u_dim_u−1,N −2





Now, the transformation to obtain a non-linear problem can really begin. As it was written in the general formulation of an optimal control problem, there were nine dierent expressions to be discretized (cost function, dynamics, initial condition, equality and inequality nal conditions, equality and inequality path constraints and two box constraints). From this point, dierent ways to discretize the dierential equations (dynamics) from the optimal control problem can be used like Euler, Middle Point, Heun, Runge-Kutta, Adams- Moulton, etc. For example, with Euler, dynamics ˙x(t) = f(t, x(t), u(t)) can be approximated as below:

x_i,j+1= x_i,j+ h. ˙x_i,j f or i = 0...dim_x− 1, f or j = 0...N − 2

which is the same as:

xi,j+1− x_i,j− h.f_i(j, x0,j, ..., xdimx−1,j, u0,j, ..., udimu−1,j) = 0, f or i = 0...dimx− 1, f or j = 0...N − 2

Whatever the method used to discretize dynamics, it gives dimx.(N − 1)dierent constraints.

In the same way, the initial condition, both inequality and equality nal conditions, may be expressed as:

xi,1− x₀(i) = 0, f or i = 0...dimx− 1 ψ1(xi,N −1) − xtf(i) ≤ 0, f or i = 0...dimx− 1 ψ₂(x_{i,N −1}) − x_t_f(i) = 0, f or i = 0...dim_x− 1

Thus, 3.dimx new constraints have just been added to the non linear problem. Moreover, the path constraints (equality and inequality) give the following constraints:

Ki(j, x0,j, ..., xdimx−1,j, u0,j, ..., udimu−1,j) ≤ 0 f or j = 0...N − 1 Ke(j, x0,j, ..., xdimx−1,j, u0,j, ..., udimu−1,j) = 0 f or j = 0...N − 1

Then, both box constraints have to be considered. Consequently, new constraints can be added:

ulower(i) ≤ u(i, j) ≤ uupper(i), f or i = 0...dimu− 1, f or j = 0...N − 2 x_lower(i) ≤ x(i, j) ≤ x_upper(i), f or i = 0...dim_x− 1, f or j = 0...N − 1

(22)

Lastly, the cost function needs to be discretized. You can observe one example of discretization below:

t_f N − 1.

N −1

X

j=0

L(x0,j, ..., xdimx−1,j, u0,j, ..., udimu−1,j) + φ(tf, xf)

Let introduce the NLP state vector z which contains both all coecients of X and U. It also includes the nal time and may be dened as:

z = (u0,0, ..., u0,N −2, ..., udimu−1,0, ..., udimu−1,N −2, ..., x0,0, ..., x0,N −1, ..., xdimx−1,0, ..., xdimx−1,N −1, tf) The dimension of the NLP state vector z is equal to dimz = (N − 1).dim_u+ N.dim_x+ 1.

The cost function can be replaced by f(z) whereas all the constraints can be replaced by g₁(z) ≤ 0and g2(z) = 0. Finally, the expression the NLP problem becomes:

min f (z) st g1(z) ≤ 0

g₂(z) = 0

This method can lead to very large non-linear problems (many variables and constraints). The way to transform an optimal control problem into a non-linear problem with the collocation will be illustrated later with a concrete example.

2.2.3.2 Direct single shooting

Let rst illustrate how the single shooting method works with the following sketch:

Figure 7: Sketch of the single shooting method

With a shooting method, only the control u is discretized. Let assume the control is piecewise constant. The objective remains to obtain a non-linear problem that might be solved by a solver. Let introduce the NLP state vector contains all coecients of U and tf. It may be dened as:

z = (u0,0, ..., u0,N −2, ..., udimu−1,0, ..., udimu−1,N −2, tf)

The dimension of the NLP state vector z is equal to dimz = (N − 1).dimu+ 1. One can notice that the dimension of z is much smaller than the state vector obtained with collocation.

This fact makes the problem more simple to solve and reduces the computation time. An integration method like Euler or Runge-Kutta is still needed to discretize dynamics. The following equations are still true with Euler:

xi,j+1− x_i,j− h.f_i(j, x0,j, ..., xdimx−1,j, u0,j, ..., udimu−1,j) = 0, f or i = 0...dim_x− 1, f or j = 0...N − 2

(23)

However, these equations are not included in the nal NLP. Starting from the initial boundary condition and using the last equation, the nal state of x composed of xN −1,0, ..., x_{N −1,dim}_x−1

may be expressed as a function of all the controls and the initial boundary. In that way, all the equations from the integration of dynamics are kept implicit. Nevertheless, the nal boundary from the optimal control problem has still to be handled. The resulted constraints in the NLP problem become:

ψ_i(xN −1,0, ..., x_{N −1,dim}_x−1) ≤ 0 inequality f inal condition ψe(xN −1,0, ..., xN −1,dimx−1) = 0 equality f inal condition

Since the nal state only depends on the controls and the initial boundary, it becomes:

ψ_i(u, x₀) ≤ 0 inequality f inal condition ψ_e(u, x₀) = 0 equality f inal condition

Both the path constraints and the cost function can be reformulated in the same way as it was made in collocation. Finally, the expression the NLP is still the same:

min f (z) st g1(z) ≤ 0

g₂(z) = 0

Compared to the collocation, the single shooting method gives smaller non-linear problems.

All the transformation with a single shooting method will be illustrated with an example later in this Master Thesis.

2.2.3.3 Direct multiple shooting

Let rst illustrate how the multiple shooting method works with the following sketch:

Figure 8: Sketch of the multiple shooting method

This method generelizes the single shooting, by adding other variables , which are guesses of the state vector x during the path and allow to do more than one shooting. Continuity constraints are added for each shooting. It leads to a greater number of constraints but gives more stability. Let nbshootings > 2be the number of shootings. The NLP vector z is composed

(24)

of controls, guesses of the state vector and tf.

z =(u_0,0, ..., u_{0,N −2}, ..., u_dim_u−1,0, ..., u_dim_u−1,N −2,

_0,0, ..., _0,nb_shootings−2, ..., _dim_x−1,0, ..., _dim_x−1,nb_shootings−2, t_f) The expression of the NLP changes but can still be expressed as:

min f (z) st g₁(z) ≤ 0

g₂(z) = 0

Three dierent methods have been explained in this part. They all allow to tranform an optimal control problem into a non-linear problem. Now, the objective is to solve the resulted NLP, which can be very dierent according the selected method and quite large in terms of variables and constraints. A powerful optimization solver is needed.

(25)

3 NUMERICAL OPTIMIZATION

3 Numerical optimization

To solve large non-linear problems, a suitable solver is mandatory. Indeed, low-thrust problems might result to NLP with thousands of variables and constraints. First, dierent state-of- the-art solvers are presented. They might be written in dierent languages. Then, a selection is performed based on dierent criteria and a test case in order to choose the most appropriate solver to solve a low-thrust problem solved by a direct method.

3.1 State of the art

The objective of this part is to list the dierent solvers which might be able to solve a large non-linear problem, as you can see in the table below.

Dierent parameters may be distinguished for a given solver like the language (modeling language: AMPL, GAMS, LING/LINGO, MPL, AIMMS, APMonitor, Mathematica, AS- CEND.. or standard language: Fortran, C, C++, Java, Python, Matlab...), the algorithm (local, global) or the licence.

Solver Algorithms (local) Algorithms(global) Language License Characteristics

Algencan AugmentedLagrangian type x Fortran.

Interfaces with C/C++

GNUGeneral Public License

Minimizing a smooth function with a

potentially large number of variables and

box-constraints COBYLA

Constrained Optimization by Linear

Approximations x

Fortran.

Interfaces with Matlab, Python, C

GNUGeneral Public License

Single precision only, minimizing a nonlinear function subject to nonlinear inequality constraints only ALGLIB

Augmented Lagrangian, AGS (Adaptive Gradient Sampling)

x C, C++ ALGLIB

Free Edition with GPL

Nonlinearly equality/inequality constrained optimization Ipopt Interior point

method x C++via a

Java interface

Eclipse Public License

Good and ecient local optimization technique

Joptimizer

Barrier

interior-point and Primal-dual interior-point methods

x Java Apache

License

Solve a convex minimization with equality and inequality constraints

NLopt Nelder-Mead

(NM), SQP various simple (21 algos)

C, C++, Fortran, Matlab

free/open- source

Good toolbox that provides many useful local optimization techniques, but less useful global techniques Artelys

Knitro

Two interior point algorithms, two active set algorithms

x C, C++,

Fortran, Java,

or Python Proprietary Nonlinear programming problems

Table 2: NLP Solvers (part 1/3)

(26)

SNOPT

Single technique (sequential quadratic programming)

x

Fortran, but interfaces with C, C++, Python and MATLAB

Proprietary Very ecient SQP algorithm

MOSEK Interior-point

optimizer x

Interfaces with C,Java and Python languages

Proprietary Convex nonlinear optimization problems

PAGMO Nelder-Mead (NM), SQP

Dierential evolution (DE), Particle swarm optimization (PSO), ASA, generic algorithms (GA)

C++ Available

Good well-structured toolbox. Well

thought-out architecture.

Includes many techniques.

WORHP Interior PointMethod x

Fortran and C implementation. Can be used from C/C++

Proprietary Continuous large scale nonlinear optimization

MIDACO x Derivative-free,

evolutionary algorithm

Java, Matlab, Python, C/C++, R and Fortran

Limited

NLP, suitable for problems with up to thousands of variables and constraints CONOPT Generalized

Reduced Gradient

(GRG) x Fortran

Available as a Fortran Subroutine library.

Large-Scale Nonlinear Optimization (NLP)

MINOS

Linearly Constrained Lagrangian Method

x Fortran Proprietary No xed limit on

problem size Lancelot Augmented

Lagrangian

Function x Fortran Academic

use Large-scale Nonlinear Optimization

BARON x Branch-and-

bound optimization

matlab + Developer- oriented parser interface

Academic andcommercial BARON licenses

Nonconvex optimization problems to global optimality

Filter SD

Linear

Complementarity Problem (LCP) Solver or Quadratic Programming (QP) Solver

x Fortran Eclipse

Public License 1.0

Smooth, twice

dierentiable, nonlinear programs.

Matlab SQP x Matlab Proprietary

Optimization function fmincon which implements a SQP scheme

GSL Nelder-Mead,

various others SA, MC C++ Available Scientic library Table 3: NLP Solvers (part 2/3)

(27)

3.2 Solver selection 3 NUMERICAL OPTIMIZATION

METSlib Single Method SA, few

others C++ Open-source

GPL

Non-mature toolbox, few and simple techniques only

OPT++ Various x C++ Public

License Provides various local optimization methods

OptSolve++Various x C++ Available Various interesting local

techniques

Tfmin

Shooting method (for the numerical solution of continuous 3D minimum time orbit transfer around the Earth)

x

Fortran and Matlab package (use the libraries Minpack and Adifor)

OpenSource To solve the TPBVP

Table 4: NLP solvers (part 3/3) 3.2 Solver selection

3.2.1 Description of the method of selection

In the previous part, a table has showed all the solvers which can be used in applied mathematics. The method of selection is based on dierent points. First, it is quite convenient to use a solver, which has an open-source code. It allows to understand the code and make some changes to better answer to a specic problem. Moreover, at Thales Services, both C++

and Java are often used. That is why the priority is given to these languages. Then, the solver should be able to handle very large non-linear problems with potentially thousands of variables and constraints. It should be able to face both equality and inequality constraints.

3.2.2 Choice of the solver

Two solvers have nally been selected: Ipopt and Nlopt. They are both open-source, coded in C++, able to face very large scale NLP with both equality and inequality constraints.

3.2.2.1 Nlopt

This solver is callable in C, C++, Fortran or Matlab. It includes dierent open-source algorithms (local or global) which often support problems with millions of variables and thousands of constraints. Some algorithms do not need to know the gradients whereas others are not able to run if the user do not provide the gradients. However, some algorithms do not support both non-linear equality and inequality constraints.

Among all of the algorithms available in Nlopt, only Cobyla (local derivative-free algorithm) and SLSQP (local gradient-based algorithm) have the characteristics to handle a non-linear problem provided by a low-thrust optimization. Cobyla is based on a constrained optimization by linear approximations. The main advantage of this algorithm is the possibility for the user to not know the gradients. On the contrary, SLSQP is based on a sequential quadratic programming. It is possible to have a better precision but you always need to compute the gradients before running this algorithm.

3.2.2.2 Ipopt

Ipopt is an open-free software package coded in C++ for large-scale linear or non-linear problems. Ipopt's algorithm is based on a interior point optimizer. It exists a Java interface

(28)

(Jipopt) which will be primordial in this report. Indeed, most of the libraries related to orbital maneuvers at Thales services are coded in Java. Jipopt allows to run Ipopt by providing a NLP in Java code. Since an interior point optimizer is used, gradients are always needed to run Ipopt. Another advantage of this solver is a great number of available options such as the tolerance on the dual infeasibility or on the constraint violation, the possibility to scale the problem etc.

3.2.3 NLP Test case

To validate the choice made in the previous part, dierent NLP problems have been tested with both solvers Nlopt and Ipopt. In this part, a rst test case is solved with Nlopt and a second one with Ipopt.

3.2.3.1 Resolution with Nlopt A very simple NLP is presented below:

min 100.(x1− x²₀)²+ (1 − x0)² st − 1.5 = x₁

This problem has only two variables and a lower boundary condition on one variable. One avantage of this test is that the analytical solution is known and can easily be compared with the numerical solution given by the solver. With an initial guess for x equal to xinitial = (−2, 1), the analytical solution is an objective function equal to 0 and an optimal point x∗

equal to (1, 1).

After computation with Nlopt (SLSQP algorithm with a tolerance equal to 10⁻⁵), the following results are found:

f (x) = 0.00173396 x∗ = (1.03482, 1.06857)

These results are close enough to validate the choice of Nlopt for a small non-linear problem.

The C++ code is available in the appendix A. As you can see in the appendix, the gradients are provided to help the SLSQP algorithm, which is a gradient-based algorithm.

Note that this NLP have similarly been solved with Ipopt during the internhip and the same results were obtained.

3.2.3.1 Resolution with Ipopt

A much more sophisticated problem is presented below:

min x1.x4(x1+ x2+ x3) + x3

st x1.x2.x3.x4≥ 25 x²₁+ x²₂+ x²₃+ x²₄ = 40 1 ≤ x1, , x2, , x3, , x4 ≤ 5

This problem has four variables, two constraints and boundary conditions for each variable.

One avantage of this test is that the analytical solution is known and can easily be compared with the numerical solution given by the solver. With an initial guess for x equal to xinitial = (1, 5, 5, 1), the minimized objective function is 0 and the optimal point x∗ is

(29)

3.3 Validation of the choice 3 NUMERICAL OPTIMIZATION

equal to (1.000, 4.7430, 3.8211, 1.3794). After computation with Ipopt, the following results are found:

f (x) = 17.014

x∗ = (1, 4.743, 3.82115, 1.37941)

These results are close enough to validate the choice of Ipopt for a more complex non-linear problem.

Note that this NLP have similarly been solved with Nlopt during the internhip and the same results were obtained.

3.3 Validation of the choice

In the previous part, both Nlopt and Ipopt eciencies have been showed on NLP test cases.

Nevertheless, in this Master Thesis, the nal objective is to handle low-thrust orbital maneuvers. Such a maneuver can be expressed thanks to an optimal control problem. In this part, an optimal control problem has been chosen: this is the Brachistochrone problem.

3.3.1 The Brachistochrone problem

The Brachistochrone problems is a famous optimal control problem which can be written as below:

min t_f

st x(t) = v(t).sin(u(t))˙

˙

y(t) = v(t).cos(u(t))

˙v(t) = g0.cos(u(t)) x(0) = 0

y(0) = 0 v(0) = 0 x(t_f) = 2 y(tf) = 2 v(tf) f ree

As you can see above, the Brachistochrone problem is a time-optimal problem with three variables: two variables for the position (x, y) and one variable for the speed v. The three

rst constraints correspond to dynamics. Then, there are both initial and nal conditions.

The aim is to nd which control u minimizes the nal time with respect to all the constraints.

This optimal control problem is interesting for this report because this is a trajectory optimization problem. In this part, the Brachistochrone problem is transformed in a non-linear problem with a direct method and nally solved with both Nlopt and Ipopt solvers. That illustrates the fact that it is possible to handle a physical trajectory optimization problem using a direct method and a suitable optimization solver. Once again, the solution of this optimal control problem is known and can be found in the scientic litterature. A comparison with the theoretical results is also made in this part. Even if the Brachistochrone problem remains easier to solve than a low-thrust problem, it a signicant rst step.

(30)

3.3.2 Transformation into NLP problem 3.3.1.1 Transformation with collocation

In this part, the Brachistochrone problem is transformed into a non-linear problem with a collocation method. The following parameters are rst introduced:

• N = number of points

• Step length: h = _{N −1}^t^f .

Dierent methods can be selected to integrate dynamics such as Euler, Runge-Kutte, Heun- Method... All of them lead to a NLP. To illustrate as simply as possible the transformation, Euler's method has been chosen in this part. Euler's method implies the following equations:

x((k + 1).h) = x(k.h) + tf

N − 1.v(k.h).sin(u(k.h)), k = 0..N − 2 y((k + 1).h) = y(k.h) + tf

N − 1.v(k.h).cos(u(k.h)), k = 0..N − 2 v((k + 1).h) = v(k.h) + tf

N − 1.g0.cos(u(k.h)), k = 0..N − 2

The three continuous constraints from dynamics have just been converted into 3.(N − 1) discrete constraints. Let uvect = [u(0), ..., u(N − 2)], xvect = [x(0), ..., x(N − 1)], yvect = [y(0), ..., y(N − 1)], vvect= [v(0), ..., v(N − 1)].

Our new problem becomes:

min tf

st xvect(0) = 0

x_vect(i + 1) − x_vect(i) − tf

N − 1.v_vect(i).sin(u_vect(i)) = 0, i = 0...N − 2 yvect(0) = 0

y_vect(i + 1) − y_vect(i) − t_f

N − 1.v_vect(i).cos(u_vect(i)) = 0, i = 0...N − 2 v_vect(0) = 0

vvect(i + 1) − vvect(i) − t_f

N − 1.g0.cos(uvect(i)) = 0, i = 0...N − 2 x_vect(N − 1) − 2 = 0

y_vect(N − 1) − 2 = 0

Let introduce the NLP variable z:

z = [u_vect, x_vect, y_vect, v_vect, t_f]

z contains all the variables from the initial problem (x, y, v) at each step, which gives a number of variables equal to: 3.N + N − 1 + 1 = 4.N for the NLP problem. The number of constraints can be computed as well: 3.(N − 1) + 3 + 2 = 3.N + 2. Finally, the NLP can be expressed as:

(31)

min z(4N − 1) st z(N − 1) = 0

z(N + i) − z(N − 1 + i) −z(4N − 1)

N − 1 .z(3N − 1 + i).sin(z(i)) = 0, i = 0..N − 2 z(2N − 1) = 0

z(2N + i) − z(2N − 1 + i) − z(4N − 1)

N − 1 .z(3N − 1 + i).cos(z(i)) = 0, i = 0..N − 2 z(3N − 1) = 0

z(3N + i) − z(3N − 1 + i) − z(4N − 1)

N − 1 .g0.cos(z(i)) = 0, i = 0..N − 2 z(2N − 2) − 2 = 0

z(3N − 2) − 2 = 0

which corresponds to a more generelized form:

min f (z) st g(z) = 0

3.3.1.2 Transformation with single shooting

In this part, only u is discretized such as u becomes piecewise constant. The same parameters N and the step length are introduced and dynamics give the same 3.(N − 1) equations as the previous part (still with Euler):

x(i + 1) − x(i) − t_f

N − 1.v(i).sin(u(i)) = 0, i = 0..N − 2 y(i + 1) − y(i) − t_f

N − 1.v(i).cos(u(i)) = 0, i = 0..N − 2 v(i + 1) − v(i) − t_f

N − 1.g₀.cos(u(i)) = 0, i = 0..N − 2

Now, a loop can be used in order to obtain x(N-1), y(N-1) and v(N-1) as functions of all the controls u(0), ... , u(N-2). Indeed, from the initial conditions x(0) = 0, y(0) = 0, v(0) = 0, the rst iteration implies:

x(1) = 0 y(1) = 0 v(1) = t_f

N − 1.g0.cos(u(0))

The last iteration allows us to express x(N-1), y(N-1) and v(N-1) as functions of all the controls u(0), ... , u(N-2) and the nal time tf. The analytic form of these functions is not given because this is a quite complex equation which does not bring more information. Let x(N − 1) = g1(u(0), ..., u(N − 2)) and y(N − 1) = g2(u(0), ..., u(N − 2)). x(N-1) and y(N-1) are constrained by the nal condition and must be qual to 2. Consequently, our new problem becomes:

min t_f

st g1(u(0), ..., u(N − 2), tf) = 2 g₂(u(0), ..., u(N − 2), t_f) = 2

(32)

The NLP state vector z can now be introduced:

z = (u(0), ..., u(N − 2), tf)

Finally, the following NLP is obtained:

min z(N − 1) st g₁(z) − 2 = 0

g₂(z) − 2 = 0

Compared to the NLP from the collocation, this problem has only two constraints and N variables.

3.3.1.3 Transformation with multiple shooting

Here, discontinuous state trajectories are allowed at intermediate NLP iterates. Extra de- cision ^k that play the role of intermediary conditions for the state variables x(t), y(t) and v(t). In this part, to simplify the equations, the tranformation is made with only three shootings. That means that there are only two intermediary points ¹ and ². It is convenient to choose a total number of points which is a multiple of the number of shootings. Let us remind that Euler's method is still used to propagate the dierential equations.

First shooting (from the boundary condition):

x(0) = 0 y(0) = 0 v(0) = 0

x(i + 1) = x(i) + t_f

N − 1.v(i).sin(u(i)), i = 0...N/3 − 1 y(i + 1) = y(i) + t_f

N − 1.v(i).cos(u(i)), i = 0...N/3 − 1 v(i + 1) = v(i) + t_f

N − 1.g₀.cos(u(i)), i = 0...N/3 − 1 Second shooting (from the rst intermediary point ¹):

x(N/3) = ¹₁ y(N/3) = ¹₂ v(N/3) = ¹₃

x(i + 1) = x(i) + t_f

N − 1.v(i).sin(u(i)), i = N/3...2N/3 − 1 y(i + 1) = y(i) + t_f

N − 1.v(i).cos(u(i)), i = N/3...2N/3 − 1 v(i + 1) = v(i) + t_f

N − 1.g₀.cos(u(i)), i = N/3...2N/3 − 1

(33)

Third shooting (from the second intermediary point ²):

x(2N/3) = ²₁ y(2N/3) = ²₂ v(2N/3) = ²₃

x(i + 1) = x(i) + t_f

N − 1.v(i).sin(u(i)), i = 2N/3...N − 2 y(i + 1) = y(i) + t_f

N − 1.v(i).cos(u(i)), i = 2N/3...N − 2 v(i + 1) = v(i) + t_f

N − 1.g₀.cos(u(i)), i = 2N/3...N − 2 x(N − 1) = 2

y(N − 1) = 2

Similarly to the single shooting method, with a loop, x(N − 1), y(N − 1), x(N/3), y(N/3), v(N/3), x(2N/3), y(2N/3), v(2N/3) can be expressed as functions of u(0), ..., u(N − 2) and tf, which are respectively noted g1, g2, g3, g4, g5, g6, g7, g8. Our problem becomes:

min tf

st x(N/3) = g₁(u(0), ..., u(N/3 − 1), t_f) = ¹₁ y(N/3) = g2(u(0), ..., u(N/3 − 1), tf) = ¹₂ v(N/3) = g₃(u(0), ..., u(N/3 − 1), t_f) = ¹₃ x(2N/3) = g4(u(N/3), ..., u(2N/3 − 1), t_f) = ²₁ y(2N/3) = g₅(u(N/3), ..., u(2N/3 − 1), t_f) = ²₂ v(2N/3) = g6(u(N/3), ..., u(2N/3 − 1), t_f) = ²₃ x(N − 1) = g7(u(2N/3), ..., u(N − 2), tf) = 2 y(N − 1) = g7(u(2N/3), ..., u(N − 2), tf) = 2

Let introduce z = (u(0), ..., u(N − 2), ¹₁, ¹₂, ¹₃, ²₁, ²₂, ²₃, tf). The dimension of z is equal to N + 6(ie N + 6 variables in the NLP). The number of constraints is equal to 8. Finally, the NLP problem can be written as:

min z(N + 5) st g₁(z) − 2 = 0

g2(z) − 2 = 0 g₃(z) − ¹₁ = 0 g₄(z) − ¹₂ = 0 g5(z) − ¹₃ = 0 g₆(z) − ²₁ = 0 g7(z) − ²₂ = 0 g₈(z) − ²₃ = 0

This problem was obtained with only three shootings. Let nbshootings be the number of shootings. More generally, if the Brachistochrone problem is transformed into a NLP with a multiple shooting, the number of variables is equal to N + 3.(nbshootings − 1) whereas the number of constraints is 2 + 3.(nbshootings− 1).

(34)

3.3.3 Final resolution

In the last part, the Brachistochrone problem has been transformed into a NLP problem, using three dierent methods (collocation, single shooting, multiple shooting). In this section, both resolutions with Nlopt and Ipopt are detailed.

3.3.3.1 Resolution with Nlopt

Here, the Brachistochrone problem is solved with Nlopt thanks to a C++ code, which implements all the functions and equations from the previous part and a nite dierence to compute the gradients. The conditions of the resolution are the following ones:

1. Collocation used to transform the optimal control problem

• Number of points N = 10

• Euler's method to integrate dynamics of the Brachistochrone problem

• Time-step for Euler's integration equal to _{N −1}^t^f . 2. SLSQP algorithm to solve the NLP problem

• Tolerance of the solver set to 10⁻⁴.

• Initial guess to start the computation:

z = (u(0), ..., u(N − 2), x(0), ..., x(N − 1), y(0), ..., y(N − 1), v(0), ..., v(N − 1), t_f) = (0, ..., 0, 0.8).

• Gradients provided with nite dierence (central dierence)

∂

∂z0

f (z0, ..., zN −1) = f (z0+ h, z1..., zN −1) − f (z0− h, z₁, ..., zN −1) 2.h

The optimal objective function (nal time) is found equal to 0.860849. The following gure illustate the optimal solution (control u) and the evolution of the main parameters (x, y, v) of the Brachistochrone problem:

Figure 9: Brachistochrone solution with Nlopt (Collocation)

(35)

These results are very satisfactory. The nal time (tf = 0.86) is slightly above the theoretical time 0.82. The control u is linear. Both evolutions of the position and the speed are very close to the expected ones from the scientic litterature. More precision can be achieved by increasing the number of points N.

3.3.3.2 Resolution with Ipopt

Here, the Brachistochrone problem is solved with Ipopt thanks to a C++ code, which implements all the functions and equations from the previous part and a nite dierence to compute the gradients. The conditions of the resolution are the following ones:

1. Single shooting used to tranform the optimal control problem

• Number of points N = 10

• Euler's method to integrate the dynamics of the Brachistochrone problem

• Time-step for Euler's integration equal to _{N −1}^t^f . 2. Interior point optimizer to solve the nlp problem

• Tolerance of the solver set to 10⁻⁴.

• Initial guess to start the computation: z = (u(0), ..., u(N − 2), tf) = (0, ..., 0, 0.8).

• Gradients provided with nite dierence (backward dierence)

∂

∂z₀f (z₀, ..., z_{N −1}) = f (z₀, z₁, ..., z_{N −1}) − f (z₀− h, z₁, ..., z_{N −1})

h (1)

The following gure illustates the optimal solution of the Brachistochrone problem with the previous conditions:

Figure 10: Brachistochrone solution with Ipopt (Single shooting)

A nal time of 0.8249 is found and the control u is still almost linear. All the curves are very similar to the ones in the last paragraph, even if the direct method (single shooting instead

Study and Industrialization of Optimization Methods for Low-Thrust Orbital Maneuvers

Study and Industrialization of Optimization Methods for Low- Thrust Orbital Maneuvers

ROMAIN BOYER

Study and Industrialization of Optimization Methods for Low-Thrust Orbital Maneuvers

Acknowledgements

Contents

List of Figures

List of Tables

1 Introduction

2 Modelling and optimization techniques

3 Numerical optimization