Optimal space launcher trajectories in case of a singular arc eric - - PowerPoint PPT Presentation

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Optimal space launcher trajectories in case of a singular arc

Fr´ ed´ eric Bonnans1, Pierre Martinon1 and Emmanuel Tr´ elat2 CEA-EDF-INRIA School Optimal Control : Algorithms and Applications May 30 - June 1st 2007

1INRIA FUTURS, team COMMANDS, and CMAP, Ecole Polytechnique 2Universit´

e d’Orl´ eans, MAPMO

1 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

2 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Problem statement - Shooting method

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

3 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Problem statement

Context: generalized Goddard problem We consider optimal 3D trajectories for a launch vehicle. Objective: payload maximization ie reach a prescribed position with a maximal final mass. Due to aerodynamic forces, possible presence of singular arcs, which means here that optimal trajectories may include arcs with a non maximal thrust.

4 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

General formulation and resolution methods

Optimal control problem (P)            Min g0(tf , x(tf )) + tf

t0 l(t, x, u)dt

Objective ˙ x = f (t, x, u) Dynamics u ∈ U Admissible controls ψ0(t0, x(t0)) = 0 Initial conditions ψ1(tf , x(tf )) = 0 Terminal conditions

• Direct methods
• state-control discretization

→ nonlinear optimization problem (NLP)

• Indirect methods
• optimality necessary conditions (Pontryagin’s Minimum Principle)

→ non linear system (shooting methods)

5 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Pontryagin’s Minimum Principle We introduce the costate p, and define the Hamiltonian: H(t, x, p, u) = l+ < p|f > Under the assumptions:

• ∃(x, u) admissible for (P), with x AC and u measurable
• f , l are C 0 wrt u, and C 1 wrt t, x
• ψ0, ψ1 are C 1 wrt x

Then any optimal pair (¯ x, ¯ u) satisfies:

• ∃¯

p = 0, AC, such that ˙ ¯ x = ∂H

∂p (t, ¯

x, ¯ p, ¯ u) ˙ ¯ p = −∂H

∂x (t, ¯

x, ¯ p, ¯ u)

• ¯

u minimizes H a.e. over [t0, tf ].

• ¯

x, ¯ p satisfy the “transversality conditions” at t0 and tf

6 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Boundary Value Problem (BVP) By noting y = (x, p), we obtain the following problem (BVP) ˙ y(t) = ϕ(t, y(t)) Boundary conditions at t0 and tf Initial value problem (IVP) By noting z the part of (x(t0), p(t0)) not set by the boundary conditions, we define (IVP) ˙ y(t) = ϕ(t, y(t)) y(t0) = y0(z) Shooting function We note ˜ y(t, z) the solution of (IVP) and define the shooting function S : z → value of boundary conditions for ˜ y(tf , z) Then solving (BVP) is equivalent to finding a zero of S.

7 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Shooting method overview Optimal control problem (P) ↓ PMP, Boundary Value Problem (BVP) Hamiltonian system, H-minimal control ↓ Initial value problem (IVP) Shooting function S ↓ Shooting method: solve S(z) = 0 (typically with a Newton method)

• Fast and accurate method
• Can be extremely sensitive to the initial point

8 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Minimum principle - Singular arcs

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

9 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Problem modelization

State (Position - Speed - Mass)      ˙ r = v ˙ v = −D(r,v)

m v v − g(r) + C u m

˙ m = −bu with (r(t), v(t), m(t)) ∈ R3 × R3 × R, D(r, v) > 0 the drag component, g(r) the gravity, C > 0 the maximal thrust and b > 0. Control u(t) ∈ R3 is the normalized control corresponding to the thrust, submitted to the constraint u(t) ≤ 1.

10 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Objective: payload maximization We maximize the final mass m(tf ), which is equivalent to minimizing the cost tf u(t) dt Initial and final conditions    (r(t0), v(t0), m(t0)) = (r0, v0, m0) r(tf ) = rf , vf is free, mf is free. Final time tf is free.

11 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Application of the Minimum principle

Hamiltonian H = (1 − bpm)u + pr, v +

• pv, −D(r, v)

m v v − g(r) + C u m

• Costate dynamics

       ˙ pr = 1

m pv,v v ∂D ∂r +

• pv, ∂g

∂r

• ˙

pv = −pr + 1

m pv,v v ∂D ∂v + D m pv v − D mpv, v v v3

˙ pm = 1

m

• pv, −D(r,v)

m v v + C u m

• 12 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Optimal control - Study of singular arcs

Optimal control u∗ minimizes H, and therefore (1 − bpm)u + C

mpv, u.

For a nondegenerate extremal:

1 The set T := {t ∈ [0, tf ] | pv(t) = 0} has a finite cardinal. 2 There exists a measurable function α on [0, tf ], with values in

[0, 1], such that u(t) = −α(t) pv(t) pv(t), a.e. on [0, tf ].

3 Define the switching function ψ(t) = 1 − bpm(t) − Cpv(t)

m(t)

. Then if ψ(t) < 0 then α(t) = 1, if ψ(t) > 0 then α(t) = 0.

13 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Singular arcs A singular arc occurs if ψ vanishes on a subset of positive measure

• f [0, tf ].

The optimal control (more precisely α(t)) is then obtained by derivating ψ = 0 until u appears explicitely.

• First derivative: ˙

ψ = 0 → u disappears (as expected)

• Second derivative: ¨

ψ = 0 → symbolic calculus yields A(r, v, m, pr, pv, pm) α = B(r, v, m, pr, pv, pm) In the numerical simulations, we check that over a singular arc

• A is non-zero so the relation is not trivial
• B

A ∈ [0, 1] so that u(t) ≤ 1 is satisfied

• A < 0 (generalized Legendre-Clebsch optimality condition).

14 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Control structure Knowing the expression of the singular control is not enough. We also need to know the control structure:

• at time t, are we on a singular arc or not ?

How to determine the number and location of singular arcs ? We regularize the problem, and try to detect the structure as the regularization tends to 0. → Continuation approach

15 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Continuation approach - Simplicial homotopy

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

16 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Continuation approach

We consider a family of problems (Pλ)λ∈[0,1], such that

• (P1) is the target problem (P)
• we known how to solve (P0).

The homotopy H connects (P0) to (P1): H : (λ, z) → Sλ(z) Continuation: follow the zero path of H, from λ = 0 to λ = 1.

• start from the known zero of H(·, 0) = S0.
• if we reach λ = 1 we have a zero of H(·, 1) = S1 = S.

Note: the existence of such a path is not automatic...

17 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Main homotopy: regularization We regularize the problem by adding a quadratic term to the criterion: Min tf u(t) + (1 − λ)u(t)2 dt, λ ∈ [0, 1] For λ < 1 the regularized Hamiltonian is strictly convex wrt u. The new control law is      if ψ(t) < −2(1 − λ) then u∗ = − pv(t)

pv(t)

if ψ(t) > 0 then u∗ = 0 else u∗ = − pv

pv pvC/m−(1−bpm) 2(1−λ)

→ no commutations or singular arc anymore

18 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Path following

• Discrete homotopy

Choose a sequence (λk) from λ = 0 to λ = 1 and solve the (Pλk) → finding a suitable sequence can be quite difficult in practice.

• Differential homotopy (Predictor-Corrector method)

Follow the zero path as a differentiable curve → fast and accurate, but requires some regularity, and uses H′

• Simplicial homotopy (Piecewise Linear method)

Builds a PL approximation of the zero path → slower and less accurate, but more robust (no derivatives)

19 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Recall the aim of the continuation approach

• predict the control structure of the solution

→ number and location of singular arcs

• provide an approximate solution

→ initialization for the shooting method with singular arcs Method choice

• high precision is not required
• numerical difficulties are expected as the regularized problem

becomes nearly singular when λ → 1. → we use a simplicial method.

20 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Simplicial method The space is triangulated into simplices, and H is interpolated over the vertices of the simplices. A simplicial algorithm builds a sequence of transverse simplices that follow the zero path of the homotopy. A transverse simplex has exactly two completely labeled faces (roughly, that contain a zero of the PL approximation of H).

x*

x0 followed zero path zero of homotopy PL approximation completely labeled face transverse simplex

Initialization: provide the first transverse simplex. Afterwards, the next transverse simplices are determined iteratively by a lexicographic minimization step and pivoting rules.

21 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Examples in R3 for two triangulations (fixed and refining)

−1.6 −1.55 −1.5 −1.45 −1.4 −1.35 −1.5 −1.4 −1.3 0.2 0.4 0.6 0.8 1 Z1(0) Z2(0) Lambda Lambda = 1 X = −1.4185 −1.4173 Simplx: 30 −1.55 −1.5 −1.45 −1.4 −1.55 −1.5 −1.45 −1.4 0.2 0.4 0.6 0.8 1 Z1(0) Z2(0) Lambda Lambda = 1 X = −1.4142 −1.4142 Simplx: 100

Principle of an adaptive triangulation meshsize

22 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Numerical Experiments

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

23 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Problem data (cf Oberle[90]) Equations are normalized with respect to r(0), m(0), and g0. Maximal thrust modulus C = 3.5; b = 7. Gravity g(r) =

g0 r3 r, with g0 = 1.

Drag D(r, v) = KDv2e−500(r−1) with KD = 310. Initial and final conditions    r0 = (0.999949994 0.0001 0.01), v0 = (0 0 0), m0 = 1, rf = (1.01 0 0), vf is free, mf is free. tf is free.

24 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Free final time We reformulate the problem as a fixed final time on via the usual transformation t = tf s, with s ∈ [0, 1] the new time tf > 0 is an additional component of the state vector, such that

• ˙

tf = 0, and tf (0), tf (1) are free,

• the associated costate satisfies ˙

ptf (s) = −H(s, x(s), p(s)). Transversality conditions        pr(1) is free pv(1) = (0 0 0) pm(1) = 0 ptf (0) = ptf (1) = 0

25 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Preliminary homotopy: atmosphere Solving the regularized problem by single shooting is not trivial. → we introduce a preliminary homotopy on the atmosphere drag: Dθ(r, v) = θKDv2e−500(r−1), θ ∈ [0, 1], θ = 0: no aerodynamic forces θ = 1: normal aerodynamic forces NB: we set λ = 0 for this preliminary homotopy.

26 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Atmosphere homotopy

Regularized problem (λ = 0)

0.5 1 0.2 0.4 0.6 0.8 Control u1 0.5 1 −2 −1.5 −1 −0.5 0 x 10

−3

Control u2 0.5 1 −0.2 −0.15 −0.1 −0.05 Control u3 0.5 1 0.5 1 ||u|| 0.5 1 −2 −1 1 Switching function ψ θ=0 θ=1

For θ = 0, the shooting method converges from the trivial starting point z0 = (0.1 0.1 0.1 0.1 0.1 0.1 0.1), and the path following reaches θ = 1 without difficulties. → We now have a solution of the regularized problem (P0).

27 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Impact of the adaptive meshsize

• red: fixed uniform meshsize
• black: adaptive meshsize
• Strong reduction of the oscillations
• Faster path following (3 times less simplices)
• Better solution at the end of the path

28 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Regularization homotopy

Solutions for a decreasing regularization (λ = 0, 0.5, 0.8):

0.5 1 0.5 1 Control u1 0.5 1 −2 −1.5 −1 −0.5 0 x 10

−3

Control u2 0.5 1 −0.2 −0.15 −0.1 −0.05 Control u3 0.5 1 0.5 1 ||u|| 0.5 1 −2 −1 1 Switching function ψ

λ=0 λ=0.5 λ=0.8

• The evolution of u and ψ strongly suggests a singular arc

29 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Remarks:

• This path following is more difficult than for the atmosphere
• homotopy. For instance, the adaptive meshsize does not perform

better than the uniform meshsize (problems at the restarts needed at the meshsize changes).

• It becomes quite hard to continue the path following above

λ = 0.8. Trajectories with erroneous bang-bang arcs appear at the vertices of the simplices, so the accuracy of the PL approximation drops. Nonetheless, the solution obtained at λ = 0.8 already gives sufficient information to solve the problem (structure + starting point for the shooting method).

30 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Solution with singular arc

Structure: we assume the presence of one singular arc. Shooting function unknown:   

• usual IVP unknown at t0
• singular arc entry and exit times: t1, t2
• state and costate value at t1, t2 : (x1, p1), (x2, p2)

Shooting function value:   

• usual terminal and transversality conditions at tf
• switching conditions at t1, t2
• state and costate matching conditions at t1, t2

31 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Switching conditions Switching conditions indicate the entry/exit of a singular arc. The expression of the singular control enforces ¨ ψ = 0. So 2 possibilities to assure that ψ = 0 over [t1, t2] are: require ψ(t1) = ψ(t2) = 0 require ψ(t1) = ˙ ψ(t1) = 0 Matching conditions x and p must be continuous at t1, t2. (˜ x(t1), ˜ p(t1)) = (x1, p1) (˜ x(t2), ˜ p(t2)) = (x2, p2)

32 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Initialization: from the solution at λ = 0.8.

0.5 1 0.5 1 Control u1 0.5 1 −4 −3 −2 −1 0 x 10

−3

Control u2 0.5 1 −0.4 −0.3 −0.2 −0.1 Control u3 0.5 1 0.5 1 ||u|| 0.5 1 −0.5 0.5 1 Switching function ψ

The expected singular arc is clearly visible.

33 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Evolution of altitude, speed and mass of the vehicle

0.2 0.4 0.6 0.8 1 0.002 0.004 0.006 0.008 0.01 0.012 ALTITUDE (||r||−||r

0||)

TIME ||r||−||r0|| 0.2 0.4 0.6 0.8 1 0.02 0.04 0.06 0.08 0.1 0.12 0.14 SPEED (||v||) TIME ||v|| 0.2 0.4 0.6 0.8 1 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 MASS TIME m

34 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Summary of the numerical resolution All tests run on a Pentium IV 2.6GHz Nonlinear solver for the shooting method: HYBRD

• Atmosphere homotopy

Numerical integration for S: Euler 25 steps Simplices crossed: 120000 Cpu time: 48 s

• Regularization homotopy (until λ = 0.8)

Numerical integration for S: Euler 25 steps Simplices crossed: 900000 Cpu time: 350 s

• Shooting method with singular arc

Numerical integration for S: RADAU5, tol. 10−6 Shooting function norm: 5 10−4 Cpu time: 17 s

35 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Solutions with a direct method

Modelization: control-only discretization NLP solver: IPOPT (A. W¨ achter and L. T. Biegler) Formulations: free thrust level and bang-bang

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 CONTROL NORM

||U|| TIME

Direct method Bang bang solution Shooting solution

• We find a solution close to the singular arc

36 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Ongoing work: A realistic launcher model

1

Problem statement - Shooting method

2

Minimum principle - Singular arcs

3

Continuation approach - Simplicial homotopy

4

Numerical Experiments

5

Ongoing work: A realistic launcher model

37 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

A more realistic launcher model

• Stage separation (mass discontinuity)

→ additional matching conditions

• Tables for some physical data

→ interpolation, derivatives, expression of singular control ?

• Low incidence during the flight (thrust/relative speed angle)

→ cone constraint for the control ?

• Regularized problem is difficult to solve, even with no atmosphere

(many trajectories lead to crashes or overconsumption) → multiple shooting ?

38 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

An exemple of trajectory with full thrust

39 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

Mass of the different parts of the launcher

MASS EVOLUTION TIME MASS Total mass Boosters ergols Stage 1 ergols Stage 2 ergols 40 / 41

Shooting method Singular arcs Simplicial homotopy Numerical Experiments Ongoing work

A few references:

H.M. Robbins, A generalized Legendre-Clebsch condition for the singular case of

• ptimal control. IBM J. of Research and Development, 11:361–372, 1967.

A.E. Bryson and Y.C. Ho, Applied optimal control, Hemisphere Publishing, New-York, 1975.

• E. Allgower and K. Georg, Numerical Continuation Methods,

Springer-Verlag, Berlin, 1990. H.J. Oberle, Numerical computation of singular control functions in trajectory

• ptimization. Journal of Guidance, Control and Dynamics, 13:153–159, 1990.
• J. Gergaud and P. Martinon, ”An application of PL continuation methods to

singular arcs problems”, Proceedings of the 12th FGS Conf. on Optimization, Avignon, 2004, Ed. A. Seeger, Lectures Notes in Economics and Mathematical Systems, Vol. 563, pp.163-186, Springer, 2006.

• F. Bonnans, P. Martinon and E. Tr´

elat, “Singular arcs in the generalized Goddard’s Problem”, RR-6157, 2007.

41 / 41