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

add tp2

parent 1665fe8b
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
+ 18
− 408
View file @ b1d5aca8
%% 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")
```
%% Output
%% 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
```
%% Output
%% 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)
```
%% Output
Iter f(x) inf-norm Step 2-norm
------ -------------- --------------
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [NaN]
* Zero: [NaN]
* Inf-norm of residuals: 0.000000
* Iterations: 0
* Convergence: true
* |x - x'| < 1.0e-08: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 1
* Jacobian Calls (df/dx): 1
p0 (solution) = NaN
%% 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
```
%% Output
%% 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)
```
%% Output
Iter f(x) inf-norm Step 2-norm
------ -------------- --------------
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [0.1, 0.1]
* Zero: [0.1, 0.1]
* Inf-norm of residuals: 0.000000
* Iterations: 0
* Convergence: true
* |x - x'| < 1.0e-08: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 1
* Jacobian Calls (df/dx): 1
p0 (solution) =
[0.1, 0.1]
%% 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.0, text("c=0", :left, :blue, 12))
```
%% Output
%% 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$ tels que $\varepsilon^2 = 0$.
On définit une fonction
$$ F(x) = \left(x_1 + x_2-1 , x_2^2 \right) $$
avec $x = (x_1, x_2)$
1. Calculer à la main $J_F$ la Jacobienne de $F$.
2. Calculer à la main $F(x + e_1 \varepsilon )$ et $F(x + e_2 \varepsilon )$ avec $(e_1, e_2)$ les vecteurs unitaires de la base canonique, et mettre le résultat sous la forme duale. Que remarquez-vous ?
3. Coder ci-dessous $F$ et $J_F$ grâce aux nombre duaux. Calculez $F(x)$ et $J_F(x)$ pour $x = \left(1, 5 \right)$. Vérifier le résultat.
**Note.**
Voici quelques fonctions utiles :
```julia
a = Dual.(x, e1) # crée le vecteur dual x + e₁ϵ
dualpart.(a) # Récupère la partie dual du vecteur a
```
%% 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