Hampath on solving optimal control problems by indirect and path - PowerPoint PPT Presentation

Hampath – on solving optimal control problems by indirect and path following methods STRUCTURAL DYNAMICAL SYSTEMS : Computational Aspects SDS2010 O. Cots in collaboration with J.-B. Caillau and J. Gergaud June 2010 O. Cots (IMB Bourgogne) Hampath June 2010 1 / 29

Motivation ORBITE FINALE 5 LEO - GOE transfer : P 0 = 11 . 625 0 r 3 −5 Mm, e 0 = 0 . 75 et i 0 = 7 ◦ 40 ORBITE INITIALE 40 20 20 Initial mass : m 0 = 1500 kg 0 0 −20 −20 −40 −40 Low thrust : T max = 0 . 1 N r 2 r 1 5 40 20 Costs r 2 0 r 3 0 Minimization of t f −20 Maximization of the final −40 −5 −50 0 50 −50 0 50 mass r 1 r 2 Minimization of the “energy” contract O. Cots (IMB Bourgogne) Hampath June 2010 2 / 29

Equations and results � t f  0 | u | dt ⇐ ⇒ max m ( t f ) 5 min 0 r 3 −5   40  40  20  20  0 ˙ = v r 0   −20 −20  = − µ 0 −40  | r | 3 r + T max u −40 ˙ v r 2  40 r 1  2 m  20  1  m ˙ = − β T max | u | 0 ( P ) r 2 r 3 0 −20 −1  −2 −40   | u | ≤ 1 −50 0 50 −40 −20 0 20 40  r 1 r 2  1  Norm of thrust    0.5   Initial and final conditions    0  t f free 0 20 40 60 80 100 120 140 time | u | : euclidian norm F IGURE : Optimal trajectory and norm of the optimal control for T max = 10 N O. Cots (IMB Bourgogne) Hampath June 2010 3 / 29

Issues and difficulties Issue The control u is discontinuous (bang-bang solution) ⇒ (BVP) with discontinuities Idea Regularize and converge to the problem. O. Cots (IMB Bourgogne) Hampath June 2010 4 / 29

Optimal control overview 1 Differential homotopy method - hampath 2 Back to the orbital transfer 3 Conclusion 4 O. Cots (IMB Bourgogne) Hampath June 2010 5 / 29

Problems solved � t f 0 f 0 ( x ( t ) , u ( t ) , λ ) dt  g ( x ( t f ) , λ ) + Min       u ( t ) ∈ U ⊂ R m x ( t ) = f ( x ( t ) , u ( t ) , λ ) ˙ , , λ ∈ R  ( P λ )   b 0 ( x ( 0 ) , λ ) = 0 ∈ R n 0  , n 0 ≤ n    b f ( x ( t f ) , λ ) = 0 ∈ R n 1  , n 1 ≤ n Definition The Hamiltonian function is H : R n × R m × R × R − × R n − → R ( x , u , λ, p 0 , p ) H ( x , u , λ, p 0 , p ) �− → H ( x , u , λ, p 0 , p ) = p 0 f 0 ( x , u , λ )+ < p , f ( x , u , λ ) > O. Cots (IMB Bourgogne) Hampath June 2010 6 / 29

Hypothesis Normal case : p 0 � = 0 ⇒ p 0 = − 1 All functions are smooth The maximization of the Hamiltonian gives a smooth function u ( x , p , λ ) = argmax w ∈ U H ( x , w , λ, p ) Definition We call true Hamiltonian the function H r ( x , λ, p ) = H ( x , u ( x , p , λ ) , λ, p ) O. Cots (IMB Bourgogne) Hampath June 2010 7 / 29

Example � 1  0 u 2 dt min   x = − λ x + u ˙  ( P λ ) x ( 0 ) = − 1    x ( 1 ) = 0 The Hamiltonian is H ( x , u , λ, p ) = − u 2 + p ( − λ x + u ) , The maximization of the Hamiltonian gives u ( x , p , λ ) = p / 2 , the true Hamiltonian is H r ( x , λ, p ) = p 2 4 − λ xp and the Principle of the Maximum of Pontryagin (PMP) does not give conditions on p ( 0 ) nor p ( 1 ) (because x ( 0 ) and x ( 1 ) are known). O. Cots (IMB Bourgogne) Hampath June 2010 8 / 29

Shooting method ( P λ ) (PMP) x = − λ x + p  ˙ 2   p = λ p ˙  ( BVP λ ) x ( 0 ) + 1 = 0   x ( 1 ) = 0  (Shooting method) � x ( 0 ) + 1 � S ( x ( 0 ) , p ( 0 ) , λ ) = = 0 x ( 1 ) with x ( 1 ) obtained from integration (it depends on x ( 0 ) , p ( 0 ) and λ ). O. Cots (IMB Bourgogne) Hampath June 2010 9 / 29

Generally ( P λ ) (PMP)  H ( x , λ, p ) = ( ∂ H r ( x ,λ, p ) , − ∂ H r ( x ,λ, p ) p ) = � (˙ x , ˙ )  ∂ p ∂ x  ( BVP λ ) b 0 ( x ( 0 ) , p ( 0 ) , λ ) = 0 ∈ R n b f ( x ( t f ) , p ( t f ) , λ ) = 0 ∈ R n   (Shooting method) � b 0 ( x 0 , p 0 , λ ) � ( S : R 2 n + 1 → R 2 n ) S ( x 0 , p 0 , λ ) = = 0 b f ( x f , p f , λ ) with ( x f , p f ) ≡ ( x ( t f , x 0 , p 0 , λ ) , p ( t f , x 0 , p 0 , λ )) = exp t f � H ( x 0 , λ, p 0 ) O. Cots (IMB Bourgogne) Hampath June 2010 10 / 29

Optimal control overview 1 Differential homotopy method - hampath 2 Back to the orbital transfer 3 Conclusion 4 O. Cots (IMB Bourgogne) Hampath June 2010 11 / 29

Introduction to the homotopy We want to compute the zeros path of the homotopic function S : R 2 n × [ 0 , 1 ] R 2 n − → ( z 0 , λ ) �− → S ( z 0 , λ ) S is nonlinear and smooth. Assuming that 0 is a regular value for S , the solution set of S ( c ) = 0 is a smooth curve. Path following techniques : Prediction-Correction (A LLGOWER AND G EORG . [2]) Path following by SVD (D IECI AND AL . [8]) hampath for optimal control [4] O. Cots (IMB Bourgogne) Hampath June 2010 12 / 29

Tangent vector computation Assuming that c ( s ) = ( z 0 ( s ) , λ ( s )) is such as c ( 0 ) = ( z 0 0 , 0 ) 1 ( z f 0 , 1 ) S ( c ( s )) = 0 2 S ′ ( c ( s )) of full rank 3 λ c ( s ) � = 0 ˙ 4 c ( s ) ˙ c ( s ) = T ( S ′ ( c ( s ))) is determined by then ˙ S ′ ( c ( s ))˙ c ( s ) = 0 z 0 1 ( z 0 0 , 0 ) | ˙ c ( s ) | =1 2 � S ′ ( c ( s )) � is of constant sign det 3 t ˙ c ( s ) O. Cots (IMB Bourgogne) Hampath June 2010 13 / 29

Differential algorithms PC [2] : an Euler step for the prediction and a Newton method for the correction. By smooth SVD [8] : smooth update of the Jacobian decomposition factors. hampath [4] uses DOPRI 5 from E. Hairer and G. Wanner [6] [7], for the numerical integration of (without any correction) : � ˙ c ( s ) = T ( S ′ ( c ( s ))) ( IVP ) c ( 0 ) = ( z 0 0 , 0 ) Until s f such as λ ( s f ) = 1 (dense output). O. Cots (IMB Bourgogne) Hampath June 2010 14 / 29

S ′ ( z 0 , λ ) S ( z 0 , λ ) = bc ( z 0 , z ( t f , z 0 , λ ) , λ ) with bc ( z 0 , z f , λ ) = ( b 0 ( z 0 , λ ) , b f ( z f , λ )) � � ∂ bc ∂ z 0 + ∂ bc ∂ z ∂ bc ∂λ + ∂ bc ∂ z ⇒ S ′ ( z 0 , λ ) = ∂ z f ∂ z 0 ∂ z f ∂λ hampath uses T APENADE [5] (from INRIA) for the automatic differentiation : ∂ bc ∂ z 0 , ∂ bc ∂ z f and ∂ bc ∂λ computation. hampath generates the variational equations to compute ∂ z ∂ z 0 and ∂ z ∂λ . O. Cots (IMB Bourgogne) Hampath June 2010 15 / 29

S ′ ( z 0 , λ ) ∂ z ∂ z 0 ( t f , z 0 , λ ) is solution of : � Y ( t ) = ∂� ˙ H ∂ z ( z ( t , z 0 , λ ) , λ ) Y ( t ) ( VAR ) Y ( 0 ) = I 2 n with ∂� ∂ z and � H H = ( ∂ H r ∂ p , − ∂ H r ∂ x ) computed by automatic differentiation (T APENADE ). Recall : z = ( x , p ) and H r is the true Hamiltonian. ∂λ ( t f , z 0 , λ ) and ∂� Same thing for ∂ z H ∂λ . O. Cots (IMB Bourgogne) Hampath June 2010 16 / 29

ssolve and expdhvfun methods expdhvfun integrates the variational equations, with any initial � Y ( t ) = ∂� ˙ H ∂ z Y ( t ) condition. It is used to check the second order Y ( 0 ) = Y 0 optimality conditions (cf [3] et [1]). ssolve uses H YBRJ (Minpack library) as NLE solver. Besides hampath integrates the variational S λ ( z 0 ) equations so that this diagram commutes S ′ λ ( z 0 ) ( DOPRI 5 has been modified) : H YBRJ Numerical integration ( IVP ) − − − − − − − − − − − → S ( z 0 , λ ) ssolve     Derivative � Derivative � Numerical integration ∂ S ( VAR ) − − − − − − − − − − − → ∂ z 0 ( z 0 , λ ) O. Cots (IMB Bourgogne) Hampath June 2010 17 / 29

hampath use User implementation (example) bc ( z 0 , z f , λ ) ✞ Subroutine bcfun(t0,tf,n,z0,zf,hf,lpar,par,nbc,bc) !Variable declarations bc = (/ z0(1)+1 , zf(1)/) end subroutine bcfun ✝ ✆ H r ( z , λ ) ✞ Subroutine hfun(t,n,z,lpar,par,H) !Variable declarations x = z(1) p = z(2) lambda = par(1) H = p**2/4 - lambda * x * p end subroutine hfun ✝ ✆ Available M ATLAB functions after compiling ssolve to solve the problem at λ fixed ; hampath to follow the zeros path ; expdhvfun to check the optimality of the solution. O. Cots (IMB Bourgogne) Hampath June 2010 18 / 29

Global diagram bcfun hfun F ORTRAN 90 bc ( z 0 , z f , λ ) H r ( z , λ ) T APENADE � H T APENADE d � dbc H ( ∂ Hr ∂ p , − ∂ Hr ∂ x ) DOPRI 5 VAR VAR DOPRI 5 DOPRI 5 S ( z 0 , λ ) S ′ ( z 0 , λ ) T ( S ′ ( z 0 , λ )) H YBRJ DOPRI 5 ssolve hampath expdhvfun M ATLAB functions shooting method Available for use continuation method variational equations O. Cots (IMB Bourgogne) Hampath June 2010 19 / 29

Optimal control overview 1 Differential homotopy method - hampath 2 Back to the orbital transfer 3 Conclusion 4 O. Cots (IMB Bourgogne) Hampath June 2010 20 / 29

Homotopy chosen � t f | u | dt Min 0 ⇓ � t f Min ( 1 − λ ) t f + λ | u | − ( 1 − λ ) ( log | u | + log ( 1 − | u | )) dt 0 For λ = 0 the final time t f is minimized ; For λ = 1 the final mass is maximized ; For λ ∈ [ 0 , 1 ) the control u ( x , p , λ ) is smooth. O. Cots (IMB Bourgogne) Hampath June 2010 21 / 29

Zeros path - λ vs costate 1 0.5 0 −2 −1 0 1 2 3 4 5 1 0.5 0 −150 −100 −50 0 50 100 1 0.5 0 −3 −2 −1 0 1 2 3 1 0.5 0 −50 −45 −40 −35 −30 −25 −20 −15 1 0.5 0 −1 −0.5 0 0.5 1 1.5 2 2.5 1 0.5 0 −3 −2 −1 0 1 2 3 4 5 6 O. Cots (IMB Bourgogne) Hampath June 2010 22 / 29