An improved switching detection for an optimal control problem with - - PowerPoint PPT Presentation

An improved switching detection for an optimal control problem with control discontinuities

Pierre Martinon, Joseph Gergaud Conference in honour of E. Hairer’s 60th birthday Gen` eve, 17-20 June 2009

1 / 20

Orbital transfer problem Initial orbit: GTO, elliptic, small inclination Final orbit: GEO, circular, equatorial

2 / 20

Low thrust (electro-ionic) propulsion: many revolutions

−60 −40 −20 20 40 −40 −20 20 40 −5 5

−40 −20 20 40 −40 −30 −20 −10 10 20 30 40 TRAJECTORY

Minimize fuel consumption: bang-bang control, many switchings

3 / 20

Optimal control problem        Min g0(x(T)) + T

0 l(x(t), u(t))dt

Objective ˙ x = f (x, u) on [0, T] Dynamics u ∈ U Admissible controls x(0) = x0, x(T) ∈ T Initial / final conditions

• Direct methods: discretization → NLP problem

+ robust and easy to implement

• solution accuracy may be too low
• Indirect methods: Pontryagin’s Minimum Principle

+ fast and precise in good conditions

• requires a good starting point, more mathematical work

4 / 20

Define the costate p and the Hamiltonian H(x, p, u) = l(x, u)+ < p, f (x, u) > Pontryagin’s Minimum Principle: For any optimal pair (x∗, u∗), there exists p∗ = 0 such that (i) ˙ p∗ = −∂H

∂x (x∗, p∗, u∗) (adjoint equation)

(ii) u∗ minimizes H(x∗, p∗, ·) on U (iii) Transversality conditions hold at t=0 and t=T ex : x fixed → p is free; x free → p − ∇g0(x) = 0 Optimal control problem → Boundary Value Problem on (x, p) Shooting method: BVP → solve S(z) = 0

5 / 20

3D transfer problem: x = [r, v, m]; r, v, u ∈ R3; |u| ≤ 1    ˙ r = v ˙ v = −µ0 r

|r|3 + u mTmax

˙ m = −C |u| with C = Tmax/(Isp g0). Optimal control minimizing H Note g(x, p) = (1 − pm)C − |pv|/m the switching function. If g(x, p) > 0 then u∗ = 0. If g(x, p) < 0 then |u∗| = 1 and u∗ = −pv/|pv|. Autonomous problem with free final time: H(x∗(t), p∗(t) = 0 ∀t ∈ [0, tf ]. H indicates the overall quality of the integration.

6 / 20

Applying the shooting method

• Use of orbital coordinates instead of cartesian
• Continuation from an “energy” criterion

Min λ|u| + (1 − λ)|u|2 λ ∈ [0, 1] λ < 1: regularized problem, control is continuous λ = 1: original problem with control switchings Here we study only the shooting for λ = 1, initialized with the solution from λ = 0.

7 / 20

Basic approach: no switching detection → Let the variable step integrator cope with the control switchings dopri5 (E.Hairer and G.Wanner, Dormand-Prince explicit RK) Absolute and relative tolerances: 1E − 08 Shooting function Jacobian computed by finite differences

Tmax 5N 1N 0.5N 0.2N 0.1N

• bj (kg)

134.49 135.9 136.31 140.04 136.02 |S| 3.63E-06 7.48E-06 9.35E-03 2.80E-01 4.74E-01 S / Jac 92 / 6 118 / 7 393 / 26 169 / 10 184 / 10 cpu (s) 1 6 42 66 128

→ Poor solutions for thrusts below 1N.

8 / 20

The stepsize becomes very small at the switchings. The ratio of rejected steps is high (≈ 50%).

50 100 150 200 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

TIME h STEPSIZE 50 100 150 200 −1 1 2 3 4 5 x 10

−6

HAMILTONIAN TIME H

The Hamiltonian indicates some errors at the switchings. → Now we will detect the switchings during the integration.

9 / 20

Switching detection At the end of each integration step: check for a sign change of the switching function. Sign change: switching crossed

• break integration
• use dense ouptut to locate it
• perform switching
• resume integration.

This method is cheap and greatly improves the integration. Using a fixed step integrator becomes possible.

10 / 20

Switching detection: effect on the integration Total number of steps (accepted + rejected) is divided by 2. Ratio of rejected steps down from ≈ 50% to ≈ 10%.

Tmax 5N 1N 0.5N 0.2N 0.1N No detection 1304 6956 14284 35652 70723 Detection 604 2933 6087 14078 29520

50 100 150 200 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

TIME h STEPSIZE NO DETECTION DETECTION

The very small steps at the switchings are strongly reduced.

11 / 20

Effect on the shooting method (no detection / detection)

Tmax 5N 1N 0.5N 0.2N 0.1N

• bj (kg)

134.49 135.9 136.31 140.04 136.02 134.49 135.9 135.82 135.81 135.71 |S| 3.63E-06 7.48E-06 9.35E-03 2.80E-01 4.74E-01 7.72E-09 1.36E-07 4.31E-07 1.45E-06 3.05E-02 cpu (s) 1 6 42 66 128 1 2 7 18 86 switchings 22 119 244 599 1201

• Same or better solutions (objective)
• Better convergence (norm of shooting function)
• Faster convergence

However, solution for 0.1N is still poor.

12 / 20

The Hamiltonian is much better than without detection. However, there remain some small errors at the switchings.

50 100 150 200 −1 1 2 3 4 5 x 10

−6

HAMILTONIAN TIME H NO DETECTION DETECTION 40 50 60 70 80 90 100 110 5 5.5 6 6.5 7 7.5 x 10

−8

HAMILTONIAN TIME H

The switching point is only an approximation by the dense ouptut. → We try to obtain a more accurate switching point.

13 / 20

Switching correction Perform usual check for sign change of the switching function Sign change: switching crossed

• break integration
• use dense ouptut to locate it
• solve for switching time τ s.t.

g(x(τ), p(τ)) = 0 with x(τ), p(τ) integrated

• perform switching
• resume integration from τ.

The switching point now corresponds to an actual integration. Initial guess: approximate switching point from dense output. This method can preserve the symmetry of the integration.

14 / 20

Switching correction: effect on the integration The Hamiltonian is close to the basic detection. However, this time time the errors at the switchings seem gone.

50 100 150 200 −6 −4 −2 2 4 6 8 x 10

−8

HAMILTONIAN TIME H DETECTION CORRECTION 40 50 60 70 80 90 100 110 4.5 5 5.5 6 6.5 7 x 10

−8

HAMILTONIAN TIME H DETECTION CORRECTION CONTROL

Does it improve the shooting methods results ?

15 / 20

Switching correction: effect on the shooting method Shooting function norm

Tmax 5N 1N 0.5N 0.2N 0.1N No detection 3.63E-06 7.48E-06 9.35E-03 2.80E-01 4.74E-01 Detection 7.72E-09 1.36E-07 4.31E-07 1.45E-06 3.05E-02 Correction 8.47E-14 2.23E-12 3.21E-11 6.94E-10 1.89E-02

Cpu times (s)

Tmax 5N 1N 0.5N 0.2N 0.1N No detection 1 6 42 66 128 Detection 1 2 7 18 86 Correction 1 2 4 21 124

• Further improvement of the shooting function norm
• Cpu times stay close to the basic detection

Solution for 0.1N still not satisfying.

16 / 20

An alternative to finite differences for the Jacobian Consider the Initial Value Problem (IVP) ˙ y(t) = ϕ(y(t)) y(t0) = y0 Let y(·, y0) be solution of (IVP), then ∂y

∂y0 (tf , y0) is solution of

(VAR)

• ˙

Ψ(t) = ∂ϕ

∂y (y(t))Ψ(t)

Ψ(t0) = I Assume a switching from ϕI to ϕII at τ, then Ψ(τ +) = Ψ(τ −) + (ϕI(yτ) − ϕII(yτ))τ ′(y0) with τ ′(y0) obtained from g(yI(τ, y0)) = 0.

17 / 20

Alternate Jacobian: effect on the shooting Switching detection (finite differences / variational system)

Tmax 5N 1N 0.5N 0.2N 0.1N |S| 7.72E-09 1.36E-07 4.31E-07 1.45E-06 3.05E-02 9.45E-09 3.58E-07 6.90E-07 1.00E-06 8.09E-04 cpu (s) 1 2 7 18 86 1 2 5 13 66

Switching correction (finite differences / variational system)

Tmax 5N 1N 0.5N 0.2N 0.1N |S| 8.47E-14 2.23E-12 3.21E-11 6.94E-10 1.89E-02 4.32E-13 4.78E-12 5.62E-11 1.28E-10 3.69E-09 cpu (s) 1 2 4 21 124 1 2 11 63 Jac 4 3 3 6 15 2 2 1 4 11

• Improvement of convergence for 0.1N
• Faster cpu times (divided by 2 for switching correction, due to

less Jacobian evaluations for the Newton method)

18 / 20

Overall comparison of convergence and Cpu times

5N 1N 0.5N 0.2N 0.1N 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 SHOOTING METHOD CONVERGENCE THRUST |S| END − NO DETECTION END − DETECTION END − CORRECTION VAR − DETECTION VAR − CORRECTION

Switching correction + Jacobian computed by variational system → overall best convergence and cpu times.

19 / 20

Conclusions

• Improvement of the basic switching detection by using a more

accurate switching point than the approximation given by the dense ouptut.

• Combined with a better computation of the Jacobian, gives

faster and more accurate results for the shooting method, for problems with up to 1200 control switchings. Perspectives

• geometric aspect: adapt the switching correction for a symplectic

integrator, or symmetric variable step integrator ?

20 / 20