Skip to content
Snippets Groups Projects
Commit dce8de3f authored by Olivier Cots's avatar Olivier Cots
Browse files

foo

parent 13cce103
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
tp/simple-shooting-etu.ipynb
+ 10
− 1
View file @ dce8de3f
%% Cell type:markdown id:efcf7482 tags:
[<img src="https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/inp-enseeiht.jpg" alt="N7" height="100"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)
# Tir simple indirect
Le but est de résoudre par du tir simple indirect, deux problèmes simples de contrôle optimal dont le contrôle maximisant est lisse.
Le premier problème est en dimension 1 et possède des conditions aux limites simples tandis que le second est en dimension 2 et aura des conditions de transversalités sur le vecteur adjoint.
On propose un dernier exercice sur la dérivation de fonctions par les nombres duaux.
%% Cell type:code id:1e3057d1 tags:
``` julia
# Packages
using DifferentialEquations, Plots, NLsolve, ForwardDiff, DualNumbers, LinearAlgebra
```
%% Cell type:markdown id:1b19260f tags:
## Exercice 1
On considère le problème de contrôle optimal suivant :
$$
\left\{
\begin{array}{l}
\displaystyle \min\, \frac{1}{2} \int_0^{1} u^2(t) \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = -x(t) + \alpha x^2(t) + u(t), \quad u(t) \in \mathrm{R}, \quad t \in [0, 1] \text{ p.p.}, \\[1.0em]
x(0) = -1, \quad x(1) = 0.
\end{array}
\right.
$$
Pour $\alpha \in \{0, 1.5\}$, faire :
1. Afficher le flot du système hamiltonien associé et afficher 5 fronts aux temps $0$, $0.25$, $0.5$, $0.75$ et $1$ ainsi que quelques extrémales. On utilisera le TP précédent.
2. Coder la fonction de tir
$$
S(p_0) = \pi(z(t_f, x_0, p_0)) - x_f,
$$
où $t_f = 1$, $x_0=-1$, $x_f = 0$ et où $z=(x,p)$ et $\pi(z) = x$. Vous pouvez utiliser votre cours.
3. Afficher la fonction de tir pour $p_0 \in [0, 1]$.
4. Résoudre l'équation $S(p_0) = 0$ et tracer l'extrémale dans le plan de phase, cf. question 1. Vous donnerez une bonne initialisation au solveur de Newton sous-jacent.
%% Cell type:code id:47d35d5b-b659-4906-bd1a-61e369de876a tags:
``` julia
# Paramètres du problème
α = 0
t0 = 0
tf = 1
x0 = -1
x_target = 0; # xf
```
%% Cell type:code id:497d73ec-5b23-459e-8c32-ca13add12f5f tags:
``` julia
## Question 1: utiliser TP précédent
# ------------------------------------
# conseil : réutiliser plot_traj! et plot_flow! du TP1
#
# à compléter ici
# fin à compléter ici
# ------------------------------------
# plot initial et paramètres du plot
plt = plot()
ymin = -0.2
ymax = 2.5
# ------------------------------------
# plot de trajectoires
#
# à compléter ici
# fin à compléter ici
# ------------------------------------
# ------------------------------------
# plot de fronts
#
# à compléter ici
# fin à compléter ici
# ------------------------------------
# additional plots
plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax), legend = false)
plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)
plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(900, 600), xlabel = "x", ylabel = "p")
```
%% Cell type:code id:fb462dc0 tags:
``` julia
# Fonction qui construit le flot à partir d'un système hamiltonien, pour résoudre
#
# (dx/dt, dp/dt) = Hv(x, p), (x(t0), p(t0)) = (x0, p0)
#
# où Hv = (dH/dp, -dH/dx).
#
function Flow(Hv)
#
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[1:n], dz[n+1:2n] = Hv(z[1:n], z[n+1:2n])
end
# cas vectoriel
function f(tspan::Tuple{Real, Real}, x0::Vector{<:Real}, p0::Vector{<:Real}; abstol=1e-12, reltol=1e-12, saveat=0.01)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0::Real, x0::Vector{<:Real}, p0::Vector{<:Real}, tf::Real; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol(tf)[1:n], sol(tf)[n+1:2n]
end
# cas scalaire
function rhs_scalar!(dz, z, dummy, t)
dz[1], dz[2] = Hv(z[1], z[2])
end
function f(tspan::Tuple{Real, Real}, x0::Real, p0::Real; abstol=1e-12, reltol=1e-12, saveat=0.01)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs_scalar!, z0, tspan)
sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0::Real, x0::Real, p0::Real, tf::Real; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
return sol(tf)[1], sol(tf)[2]
end
return f
end;
```
%% Cell type:code id:df2c37e6 tags:
``` julia
## Question 2 : coder la fonction de tir
# Système hamiltonien
function Hv(x, p)
"""
Calcul de la dynamique hamiltonienne
-------
parameters (input)
------------------
x : état
p : état adjoint
returns
-------
ẋ, ṗ
"""
dx = 0
dp = 0
return dx, dp
end
# flot du système hamiltonien
f = Flow(Hv)
# fonction de tir
function shoot(p0)
return 0
end;
```
%% Cell type:code id:ba8b9913-f4f7-4a2f-9003-8179bf157aea tags:
``` julia
## Question 3 : afficher la fonction de tir
# Plot initial
plt_2 = plot(size=(600, 400),xlabel = "p₀", ylabel = "S(p₀)", xticks = 0:0.1:1, xlims = (0,1), ylims = (-0.5,1) )
plot!(plt_2, [0,1], [0,0], legend = false, color = :black)
p0 = range(0, 1, length = 50)
#
# compléter à partir d'ici
```
%% Cell type:code id:ea2ab41c-44d1-4e24-a99a-32f7e2a0b75a tags:
``` julia
## Question 4 : résoudre S(p0) = 0
# transformation de la fonction de tir en une fonction vectorielle
S(p0) = [shoot(p0[1])]
# jacobienne de la fonction de tir
JS(p0) = ForwardDiff.jacobian(S, p0)
# initialisation de la fonction de tir
p0_guess = [NaN] # A MODIFIER
# résolution de shoot(p0) = 0. L'itéré initial est appelé p0_guess
sol = nlsolve(S, JS, p0_guess; xtol=1e-8, method=:trust_region, show_trace=true);
println("\n", sol, "\n")
# récupération de la solution
p0_sol = sol.zero[1]
#
println("p0 (solution) = ", p0_sol)
```
%% Cell type:code id:fe60edc8-df13-454e-a12d-700ec7f2419c tags:
``` julia
# Affichage de la BC-extrémale
#
# CONSEIL : utiliser la fonction f
#
plt # A REMPLACER
```
%% Cell type:markdown id:15fe1891 tags:
## Exercice 2
On considère le problème de contrôle optimal suivant :
$$
\left\{
\begin{array}{l}
\displaystyle \min\, \frac{1}{2} \int_0^{t_f} ||u(t)||^2 \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = (0, \alpha) + u(t), \quad u(t) \in \mathrm{R}^2, \quad t \in [0, t_f] \text{ p.p.}, \\[1.0em]
x(0) = a, \quad x(t_f) \in b + \mathrm{R} v,
\end{array}
\right.
$$
avec $a = (0,0)$, $b = (1,0)$, $v = (0, 1)$ et $t_f = 1$.
Pour différente valeur de $\alpha \in [0, 2]$, faire :
1. Résoudre le problème de contrôle optimal par du tir simple.
2. Tracer la trajectoire solution dans le plan $(x_1, x_2)$.
Tracer pour la solution, attaché au point $x(t_f)$, les vecteurs $\dot{x}(t_f)$ et $p(t_f)$. Commentaires ?
**Remarque.**
Utiliser la fonction suivante pour tracer un vecteur $v = (v_1, v_2)$ aux coordonnées $(a, b)$
``` julia
quiver!(plt, [a], [b], quiver=([v1], [v2]), c="green", linewidth=2)
```
%% Cell type:code id:2a0f0b8b-5092-41cb-bb53-815c99ea95ee tags:
``` julia
# Paramètres du problème
α = 2
t0 = 0
tf = 1
x0 = [0, 0];
```
%% Cell type:code id:90248383 tags:
``` julia
## Question 1 : résolution de S(p0) = 0
# système hamiltonien
function Hv(x, p)
return [0, 0], [0, 0]
end
# flot du système hamiltonien
z = Flow(Hv)
# fonction de tir
function shoot(p0)
return [0, 0]
end
# jacobienne de la fonction de tir
jshoot(p0) = ForwardDiff.jacobian(shoot, p0)
# initialisation de la fonction de tir
p0_guess = [0.1, 0.1]
# résolution de shoot(p0) = 0. L'itéré initial est appelé p0_guess
sol = nlsolve(shoot, jshoot, p0_guess; xtol=1e-8, method=:trust_region, show_trace=true);
println("\n", sol, "\n")
# récupération de la solution
p0_sol = sol.zero
#
println("p0 (solution) = ", p0_sol)
```
%% Cell type:code id:8bd070d1-8474-4fbb-b0d5-054fdc5bffcd tags:
``` julia
## Question 2 : affichage de la solution
# calcul de la solution
ode_sol = z((t0, tf), x0, p0_sol);
# affichage de la contrainte terminale
plt = plot([1, 1], [-2, 4], c="blue", xlabel="x₁", ylabel="x₂", linewidth=2.0, xlims=(0,2.5), ylim = (-2,2), legend = false) # cible
annotate!(plt, 1.05, -1, text("c=0", :left, :blue, 12))
```
%% Cell type:markdown id:2be508c9-65a2-45a5-862d-314a973e4cca tags:
## Exercice 3
Vous venez d'utiliser le package `ForwardDiff` pour calculer la Jacobienne de la fonction de tir. Ce package utilise les nombres duaux afin de pouvoir diférencier automatiquement vos fonctions.
Les nombres duaux s'écrivent sous la forme $a + b\, \varepsilon$ avec $(a,b)\in \mathbb{R}^2$ et $\varepsilon^2 = 0$. Nous allons voir comment nous pouvons les utiliser pour calculer des dérivées.
Soit deux fonctions $f, g \colon \mathbb{R} \to \mathbb{R}$ dérivables, de dérivées respectives $f'$ et $g'$. On pose
$$
f(a + b\, \varepsilon) \coloneqq f(a) + f'(a)\, b\, \varepsilon
$$
et
$$
g(a + b\, \varepsilon) \coloneqq g(a) + g'(a)\, b\, \varepsilon.
$$
On a alors les propriétés suivantes. Posons $d = x + \varepsilon$, alors :
- $(f + g)(d) = (f+g)(x) + (f+g)'(x) \, \varepsilon$
- $(fg)(d) = (fg)(x) + (fg)'(x) \, \varepsilon$
- $(g \circ f)(d) = (g \circ f)(x) + (g \circ f)'(x) \, \varepsilon$
**Questions.**
1. Vérifier les propriétés ci-dessus.
2. Soit $f(x) = x^2$. Coder la fonction $f(x)$ et coder sa dérivée à l'aide des nombres duaux. Vérifier le résultat pour $x = 2$ par exemple.
3. Soit $g(y) = \cos(y)$. Coder la fonction $h(x) = g \circ f(x)$ et sa dérivée à l'aide des nombres duaux.
Que remarquez-vous ? Que devez-vous implémenter pour que cela fonctionne ? Vérifier.
4. Bonus. Soit $f(x, y) = x^2 + \cos(y)$. Coder la fonction et sa dérivée, et vérifier.
**Remarques importantes.**
- Vous utiliserez les fonctions suivantes du package `DualNumbers`.
```julia
d = Dual(x, 1) # crée le vecteur dual x + ϵ
d = Dual(x, 1) # crée le nombre dual x + ϵ
realpart(d) # récupère la partie réelle de d
dualpart(d) # récupère la partie duale de d
```
ou dans le cas vectoriel (attention au `.`) :
```julia
d = Dual.([x, y], [u, v]) # crée le vecteur dual (x, y) + (u, v) ϵ
realpart.(d) # récupère la partie réelle de d : (x, y)
dualpart.(d) # récupère la partie duale de d : (u, v)
```
- En `Julia`, nous pouvons typer les fonctions pour utiliser le "polymorphisme" appelé aussi "multiple dispatch", ce qui permet de choisir le comportement d'une fonction en fonction du type des arguments. Par exemple :
```julia
julia> fun(y::Integer) = 2y
julia> fun(y::Real) = y/2
julia> fun(2)
4
julia> fun(2.0)
1.0
```
**Attention !**
La fonction $f$ étant déjà définie plus haut, définissez-la en utilisant une majuscule :
```julia
F(x) = x^2
```
%% Cell type:code id:5adb3c27-a9fa-4097-99eb-d9feda6d66d9 tags:
``` julia
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment