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

−40 −20 20 40 −40 −20 20 40 −5 5 r1 ORBITE INITIALE ORBITE FINALE r2 r3 −50 50 −40 −20 20 40 r1 r2 −50 50 −5 5 r2 r3

contract

LEO - GOE transfer : P0 = 11.625 Mm, e0 = 0.75 et i0 = 7◦ Initial mass : m0 = 1500 kg Low thrust : Tmax = 0.1 N Costs Minimization of tf Maximization of the final mass Minimization of the “energy”

• O. Cots (IMB Bourgogne)

Hampath June 2010 2 / 29

Equations and results

(P)                                min tf

0 |u|dt ⇐

⇒ max m(tf ) ˙ r = v ˙ v = − µ0

|r|3 r + Tmax u m

˙ m = −βTmax|u| |u| ≤ 1 Initial and final conditions tf free |u| : euclidian norm

−40 −20 20 40 −40 −20 20 40 −5 5 r1 r2 r3 −50 50 −40 −20 20 40 r1 r2 −40 −20 20 40 −2 −1 1 2 r2 r3 20 40 60 80 100 120 140 0.5 1 time Norm of thrust

FIGURE: Optimal trajectory and norm

• f the optimal control for Tmax = 10N
• 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

1

Optimal control overview

2

Differential homotopy method - hampath

3

Back to the orbital transfer

4

Conclusion

• O. Cots (IMB Bourgogne)

Hampath June 2010 5 / 29

Problems solved

(Pλ)                Min g(x(tf ), λ) + tf

0 f 0(x(t), u(t), λ)dt

˙ x(t) = f(x(t), u(t), λ) , u(t) ∈ U ⊂ Rm , λ ∈ R b0(x(0), λ) = 0 ∈ Rn0 , n0 ≤ n bf (x(tf ), λ) = 0 ∈ Rn1 , n1 ≤ n

Definition

The Hamiltonian function is H : Rn × Rm × R × R− × Rn − → R (x, u, λ, p0, p) − → H(x, u, λ, p0, p) H(x, u, λ, p0, p) = p0f 0(x, u, λ)+ < p, f(x, u, λ) >

• O. Cots (IMB Bourgogne)

Hampath June 2010 6 / 29

Hypothesis

Normal case : p0 = 0 ⇒ p0 = −1 All functions are smooth The maximization of the Hamiltonian gives a smooth function u(x, p, λ) = argmaxw∈UH(x, w, λ, p)

Definition

We call true Hamiltonian the function Hr(x, λ, p) = H(x, u(x, p, λ), λ, p)

• O. Cots (IMB Bourgogne)

Hampath June 2010 7 / 29

Example

(Pλ)        min 1

0 u2dt

˙ x = −λx + u x(0) = −1 x(1) = 0 The Hamiltonian is H(x, u, λ, p) = −u2 + p(−λx + u), The maximization of the Hamiltonian gives u(x, p, λ) = p/2, the true Hamiltonian is Hr(x, λ, p) = p2

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) (BVPλ)        ˙ x = −λx + p

2

˙ p = λp x(0) + 1 = 0 x(1) = 0 (Shooting method) S(x(0), p(0), λ) = x(0) + 1 x(1)

• = 0

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) (BVPλ)      (˙ x, ˙ p) = H(x, λ, p) = ( ∂Hr(x,λ,p)

∂p

, − ∂Hr(x,λ,p)

∂x

) b0(x(0), p(0), λ) = 0 ∈ Rn bf (x(tf ), p(tf ), λ) = 0 ∈ Rn (Shooting method) S(x0, p0, λ) = b0(x0, p0, λ) bf (xf , pf , λ)

• = 0

(S : R2n+1 → R2n) with (xf , pf ) ≡ (x(tf , x0, p0, λ), p(tf , x0, p0, λ)) = exptf H(x0, λ, p0)

• O. Cots (IMB Bourgogne)

Hampath June 2010 10 / 29

1

Optimal control overview

2

Differential homotopy method - hampath

3

Back to the orbital transfer

4

Conclusion

• 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 : R2n × [0, 1] − → R2n (z0, λ) − → S(z0, λ) 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 (ALLGOWER AND GEORG. [2]) Path following by SVD (DIECI AND AL. [8]) hampath for optimal control [4]

• O. Cots (IMB Bourgogne)

Hampath June 2010 12 / 29

Tangent vector computation

Assuming that c(s) = (z0(s), λ(s)) is such as

1

c(0) = (z0

0, 0)

2

S(c(s)) = 0

3

S′(c(s)) of full rank

4

˙ c(s) = 0 then ˙ c(s) = T(S′(c(s))) is determined by

1

S′(c(s))˙ c(s) = 0

2

|˙ c(s)|=1

3

det

• S′(c(s))

t˙

c(s)

• is of constant sign

λ z0 (zf

0, 1)

(z0

0, 0)

˙ 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 DOPRI5 from E. Hairer and G. Wanner [6] [7], for the numerical integration of (without any correction) : (IVP) ˙ c(s) = T(S′(c(s))) c(0) = (z0

0, 0)

Until sf such as λ(sf ) = 1 (dense output).

• O. Cots (IMB Bourgogne)

Hampath June 2010 14 / 29

S′(z0, λ)

S(z0, λ) = bc(z0, z(tf , z0, λ), λ) with bc(z0, zf , λ) = (b0(z0, λ) , bf (zf , λ)) ⇒ S′(z0, λ) =

• ∂bc

∂z0 + ∂bc ∂zf ∂z ∂z0 ∂bc ∂λ + ∂bc ∂zf ∂z ∂λ

• hampath uses TAPENADE [5] (from INRIA) for the automatic

differentiation : ∂bc

∂z0 , ∂bc ∂zf and ∂bc ∂λ computation.

hampath generates the variational equations to compute ∂z

∂z0 and ∂z ∂λ.

• O. Cots (IMB Bourgogne)

Hampath June 2010 15 / 29

S′(z0, λ)

∂z ∂z0 (tf , z0, λ) is solution of :

(VAR)

• ˙

Y(t) = ∂

H ∂z (z(t, z0, λ), λ)Y(t)

Y(0) = I2n with ∂

H ∂z and

H = ( ∂Hr

∂p , − ∂Hr ∂x ) computed by automatic

differentiation (TAPENADE). Recall : z = (x, p) and Hr is the true Hamiltonian. Same thing for ∂z

∂λ(tf , z0, λ) and ∂ H ∂λ .

• O. Cots (IMB Bourgogne)

Hampath June 2010 16 / 29

ssolve and expdhvfun methods

• ˙

Y(t) = ∂

H ∂z Y(t)

Y(0) = Y0 Sλ(z0) S′

λ(z0)

ssolve HYBRJ

expdhvfun integrates the variational equations, with any initial

• condition. It is used to check the second order
• ptimality conditions (cf [3] et [1]).

ssolve uses HYBRJ (Minpack library) as NLE solver. Besides hampath integrates the variational equations so that this diagram commutes (DOPRI5 has been modified) : (IVP)

Numerical integration

− − − − − − − − − − − → S(z0, λ)

Derivative

 

• 

 Derivative (VAR)

Numerical integration

− − − − − − − − − − − →

∂S ∂z0 (z0, λ)

• O. Cots (IMB Bourgogne)

Hampath June 2010 17 / 29

hampath use

User implementation (example) bc(z0, zf , λ)✞ Subroutine bcfun(t0,tf,n,z0,zf,hf,lpar,par,nbc,bc) !Variable declarations bc = (/ z0(1)+1 , zf(1)/) end subroutine bcfun

✝ ✆

Hr(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 MATLAB 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 bc(z0, zf , λ) hfun Hr(z, λ) dbc

• H

( ∂Hr

∂p , − ∂Hr ∂x )

S(z0, λ) d H S′(z0, λ) T(S′(z0, λ)) ssolve shooting method hampath continuation method expdhvfun variational equations TAPENADE TAPENADE

DOPRI5

VAR

DOPRI5 DOPRI5

VAR

DOPRI5

HYBRJ MATLAB functions Available for use FORTRAN 90

• O. Cots (IMB Bourgogne)

Hampath June 2010 19 / 29

1

Optimal control overview

2

Differential homotopy method - hampath

3

Back to the orbital transfer

4

Conclusion

• O. Cots (IMB Bourgogne)

Hampath June 2010 20 / 29

Homotopy chosen

Min tf |u|dt ⇓ Min (1 − λ)tf + λ tf |u| − (1 − λ) (log |u| + log(1 − |u|)) dt For λ = 0 the final time tf 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

−2 −1 1 2 3 4 5 0.5 1 −150 −100 −50 50 100 0.5 1 −3 −2 −1 1 2 3 0.5 1 −50 −45 −40 −35 −30 −25 −20 −15 0.5 1 −1 −0.5 0.5 1 1.5 2 2.5 0.5 1 −3 −2 −1 1 2 3 4 5 6 0.5 1

• O. Cots (IMB Bourgogne)

Hampath June 2010 22 / 29

The control

|u| for λ ∈ {0, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999}

−40 −20 20 40 −40 −20 20 −5 5 r1 r2 r3

λ = 0

−40 −20 20 40 −40 −20 20 −5 5 r1 r2 r3

λ = 0.9

• O. Cots (IMB Bourgogne)

Hampath June 2010 23 / 29

The control

−40 −20 20 40 −40 −20 20 40 −5 5 r1 r2 r3

λ = 0.999

• O. Cots (IMB Bourgogne)

Hampath June 2010 23 / 29

Locally optimality (assuming tf fixed)

expt,x0,λ : Rn − → Rn p0 − → x(t, x0, p0, λ) (x0, p0) is a locally optimal solution ⇔ ∀t ∈ (0, tf ] rank( d

dp0 expt,x0,λ(p0)) = n

⇔ ∀t ∈ (0, tf ]

∂x ∂p0 (t, x0, p0, λ) ∈ Rn×n is

• f full rank

10 20 30 40 50 60 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10

−3

l σ1

Smallest singular value w.r.t. L for λ ∈ {0, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999}

• O. Cots (IMB Bourgogne)

Hampath June 2010 24 / 29

1

Optimal control overview

2

Differential homotopy method - hampath

3

Back to the orbital transfer

4

Conclusion

• O. Cots (IMB Bourgogne)

Hampath June 2010 25 / 29

Points to remember

hampath is used in the case of smooth functions ; permits to treat discontinuous optimal problems by smooth regularization ; computes accurately the jacobian S′(z0, λ) ; is easy to use since it has a MATLAB interface.

• O. Cots (IMB Bourgogne)

Hampath June 2010 26 / 29

Other applications

Quantum mechanic problems The particles are controlled by laser fields and interact with the environment.

homotopy on physical parameters (characteristics of the particle). regularization by homotopy (to be smooth).

3 bodies problem The satellite goes from the Earth to the moon.

homotopy on the influence of the moon. homotopy on the maximal thrust. ...

• O. Cots (IMB Bourgogne)

Hampath June 2010 27 / 29

Improvement

Discontinuous case hampath can be applied with smooth functions and permits to have a good approximation of a discontinuous problem.

Associate the code SHOOT PACKAGE [9] (shooting method for the discontinuous case) with hampath. Study the second order optimality conditions in the discontinuous case (no complete results for now).

Continuation method and bifurcations hampath computes at each step the jacobian S′(z0, λ) used to find the tangent vector. We want to integrate the work of DIECI AND AL. [8] to have a path following by smooth SVD.

Update of the decomposition of the jacobian. Detection of bifurcation points.

• O. Cots (IMB Bourgogne)

Hampath June 2010 28 / 29

References

• A. Agrachev and Y. Sachkov.

Control theory from the geometric viewpoint, volume 87. Springer, 2004. E.L. Allgower and K. Georg. Introduction to numerical continuation methods, volume 45 of Classics In Applied Mathematics. SIAM, 2003. J.-B. Caillau B. Bonnard and E. Trélat. COTCOT http://apo.enseeiht.fr/cotcot/. J.B. Caillau, O. Cots, and J. Gergaud. HAMPATH apo.enseeiht.fr/hampath.

• A. Dervieux, L. Hascoet, and V. Pascal.

TAPENADE http://www-sop.inria.fr/tropics/tapenade.html.

• E. Hairer, S.P

. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I, Nonstiff Problems, volume 8 of Springer Serie in Computational Mathematics. Springer-Verlag, second edition, 1993.

• E. Hairer and G. Wanner.

DOPRI5 http://www.unige.ch/~hairer/prog/nonstiff/dopri5.f.

• M. G. Gasparo L. Dieci and A. Papini.

Path following by svd. Computational Science - ICCS 2006, Volume 3994/2006 :677–684, 2006. P . Martinon and J. Gergaud. SHOOT PACKAGE apo.enseeiht.fr/shoot.

• O. Cots (IMB Bourgogne)

Hampath June 2010 29 / 29