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 TRAJECTORY 40 30 5 20 0 −5 10 40 0 20 −10 0 40 −20 20 −20 0 −30 −20 −40 −40 −40 −60 −40 −20 0 20 40 Minimize fuel consumption : bang-bang control, many switchings 3 / 20

Optimal control problem � T  Min g 0 ( x ( T )) + 0 l ( x ( t ) , u ( t )) dt Objective   x = f ( x , u ) on [0 , T ] ˙ Dynamics  u ∈ U Admissible controls    x (0) = x 0 , 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 p ∗ = − ∂ H ∂ x ( x ∗ , p ∗ , u ∗ ) ( adjoint equation ) ( i ) ˙ ( 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 − ∇ g 0 ( 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 ∈ R 3 ; | u | ≤ 1  ˙ = r v  − µ 0 r | r | 3 + u v ˙ = m T max  m ˙ = − C | u | with C = T max / ( Isp g 0 ). Optimal control minimizing H Note g ( x , p ) = (1 − p m ) C − | p v | / m the switching function. If g ( x , p ) > 0 then u ∗ = 0. If g ( x , p ) < 0 then | u ∗ | = 1 and u ∗ = − p v / | p v | . Autonomous problem with free final time: H ( x ∗ ( t ) , p ∗ ( t ) = 0 ∀ t ∈ [0 , t f ]. 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: 1 E − 08 Shooting function Jacobian computed by finite differences T max 5N 1N 0.5N 0.2N 0.1N obj (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%). STEPSIZE −6 HAMILTONIAN 5 x 10 1 10 0 10 4 −1 10 3 −2 10 h H 2 −3 10 1 −4 10 −5 0 10 −6 10 −1 0 50 100 150 200 0 50 100 150 200 TIME TIME 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%. T max 5N 1N 0.5N 0.2N 0.1N No detection 1304 6956 14284 35652 70723 Detection 604 2933 6087 14078 29520 STEPSIZE 1 10 0 10 −1 10 The very small steps at the −2 10 h switchings are strongly reduced. −3 10 −4 10 −5 10 NO DETECTION DETECTION −6 10 0 50 100 150 200 TIME 11 / 20

Effect on the shooting method ( no detection / detection) T max 5N 1N 0.5N 0.2N 0.1N obj (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. HAMILTONIAN HAMILTONIAN −6 −8 5 x 10 x 10 NO DETECTION 7.5 DETECTION 4 7 3 6.5 H 2 H 6 1 5.5 0 5 −1 0 50 100 150 200 40 50 60 70 80 90 100 110 TIME TIME 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. HAMILTONIAN HAMILTONIAN −8 −8 8 x 10 x 10 6 7 4 6.5 2 6 H H 0 5.5 −2 5 DETECTION −4 DETECTION CORRECTION 4.5 CORRECTION CONTROL −6 0 50 100 150 200 40 50 60 70 80 90 100 110 TIME TIME Does it improve the shooting methods results ? 15 / 20

Switching correction: effect on the shooting method Shooting function norm T max 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) T max 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 � ˙ y ( t ) = ϕ ( y ( t )) ( IVP ) y ( t 0 ) = y 0 Let y ( · , y 0 ) be solution of ( IVP ), then ∂ y ∂ y 0 ( t f , y 0 ) is solution of � Ψ( t ) = ∂ϕ ˙ ∂ y ( y ( t ))Ψ( t ) ( VAR ) Ψ( t 0 ) = I Assume a switching from ϕ I to ϕ II at τ , then Ψ( τ + ) = Ψ( τ − ) + ( ϕ I ( y τ ) − ϕ II ( y τ )) τ ′ ( y 0 ) with τ ′ ( y 0 ) obtained from g ( y I ( τ, y 0 )) = 0. 17 / 20

Alternate Jacobian: effect on the shooting Switching detection ( finite differences / variational system) T max 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) T max 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 0 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 SHOOTING METHOD CONVERGENCE 0 10 END − NO DETECTION END − DETECTION −2 10 END − CORRECTION VAR − DETECTION −4 VAR − CORRECTION 10 −6 10 |S| −8 10 −10 10 −12 10 −14 10 5N 1N 0.5N 0.2N 0.1N THRUST 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