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

cleaning

parent 578b17ca
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
README.md
+ 0
− 10
View file @ c4af34c8
......@@ -14,16 +14,6 @@ Systèmes dynamiques
* [Sujet 1 : tp/ode.ipynb](tp/ode.ipynb)
Méthodes indirectes et directes en contrôle optimal
* [Sujet 2 : tp/simple-shooting.ipynb](tp/simple-shooting.ipynb) - Tir simple indirect
* [Sujet 3 : tp/direct-indirect.ipynb](tp/direct-indirect.ipynb) - Méthode directe et tir avec structure sur le contrôle
Calcul différentiel et équation différentielle ordinaire
* [Sujet 4 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul numérique de dérivées
* [Sujet 5 : tp/runge-kutta.ipynb](tp/runge-kutta.ipynb) - Méthodes de Runge-Kutta
**Notes de cours supplémentaires - ressources**
* [Automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf)
......
tp/derivative.ipynb deleted 100644 → 0
+ 0
− 495
View file @ 578b17ca
%% Cell type:markdown id: tags:
[<img src="https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png" alt="N7" height="100"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)
# Calcul de dérivées
- Date : 2023-2024
- Durée approximative : 1h15
## Introduction
Il existe plusieurs façon de calculer une dérivée sur un calculateur :
- par différences finies;
- par différentiation complexe;
- en utilisant la différentiation automatique;
- en utilisant le calcul formel et un générateur de code.
Nous allons étudier ici quelques cas.
On notera $\Vert{\cdot}\Vert$ la norme euclidienne usuelle et $\mathcal{N}(0,1)$ une variable aléatoire Gaussienne centrée réduite.
%% Cell type:code id: tags:
``` julia
# load packages
using DualNumbers
using DifferentialEquations
using ForwardDiff
using LinearAlgebra
using Plots
using Plots.Measures
using Printf
```
%% Cell type:markdown id: tags:
## Dérivées par différences finies avant
Soit $f$ une fonction lisse de $\mathbb{R}^{n}$ dans $\mathbb{R}^{m}$, $x$ un point de $\mathbb{R}^{n}$ et $v$ un vecteur de $\mathbb{R}^{n}$. Si on note $g \colon h \mapsto f(x+hv)$, alors on a d'après la formule de Taylor-Young :
$$
g(h) = \sum_{i=0}^{n} \frac{h^i}{i!} g^{(i)}(0) + R_n(h), \quad R_n(h) = o(h^n),
$$
ou d'après Taylor-Lagrange,
$$
\| R_n(h) \| \leq \frac{M_n h^{n+1}}{(n+1)!},
$$
de même,
$$
g(-h) = \sum_{i=0}^{n} \frac{(-h)^i}{i!} g^{(i)}(0) + R_n(h).
$$
La méthode des différences finies avants consiste à approcher la différentielle de $f$ en $x$ dans la direction $v$ par la formule suivante :
$$
\frac{f(x+hv) - f(x)}{h} =
\frac{g(h)-g(0)}{h} = g'(0) + \frac{h}{2} g^{2}(0) + \frac{h^2}{6} g^{(3)}(0) + o(h^2).
$$
L'approximation ainsi obtenue de $ g'(0) = f'(x) \cdot v \in \mathbb{R}^m$ est d'ordre 1 si $g^{(2)}(0) \neq 0$ ou au moins d'ordre 2 sinon.
**Remarque.** Sur machine, Les calculs se font en virgule flottante. On note epsilon machine, le plus petit nombre $\mathrm{eps}_\mathrm{mach}$ tel que $1+\mathrm{eps}_\mathrm{mach}\ne 1$. Cette quantité dépend des machines et de l'encodage des données. Pour l'optenir en `Julia` il suffit de taper `eps()`. On peut indiquer la précision numérique : `eps(Float64)` ou `eps(Float64(1))` pour les flottants codés sur 64 bits et `eps(Float32(1))` sur 32 bits.
Notons $\mathrm{num}(g,\, h)$ la valeur de $g(h)$ calculée numériquement et supposons que l'on puisse majorer l'erreur relative numérique par :
$$
\left\| \mathrm{num}(g,h) - g(h) \right\| := \| e_h\| \leq \mathrm{eps}_\mathrm{mach} L_f,
$$
ou $L_f$ est une constante qui dépend de la valeur de $f$ sur le domaine d'intérêt. Ainsi on a :
\begin{align*}
\left\| \frac{\mathrm{num}(g,h) - \mathrm{num}(g,0)}{h} - g'(0) \right\|
&= \left\| \frac{g(h) + e_h - g(0) - e_0}{h} - g'(0) \right\|, \\[1em]
&= \left\| \frac{R_1(h)}{h} + \frac{e_h - e_0}{h} \right\|, \\[1em]
&\leq \left\| \frac{R_1(h)}{ h} \right\| + \left\| \frac{e_h - e_0}{h} \right\|, \\[1em]
& \leq
\underbrace{ \frac{M_1 h}{2}}_{{\text{Erreur d'approximation}}} +
\underbrace{2 \frac{\mathrm{eps}_\mathrm{mach}L_f}{h}}_{{\text{Erreur numérique}}}.
\end{align*}
Le majorant trouvé atteint son minimum en
$$
h_{*} = 2 \sqrt{\frac{\mathrm{eps}_\mathrm{mach} L_f }{M_1}}.
$$
En considérant que $L_f \simeq M_1$, alors le choix se révélant le plus optimal est
$$
h_{*} \approx \sqrt{\mathrm{eps}_\mathrm{mach}}.
$$
%% Cell type:markdown id: tags:
## Dérivées par différences finies centrées
On peut utiliser un schéma de différences finies centrée pour calculer la dérviée de $g$.
$$
\frac{f(x+hv) - f(x-hv)}{2h} = \frac{g(h) - g(-h)}{2h} =
g'(0) + g^{(3)}(0) \frac{h^2}{6} + \mathcal{O}(h^4),
$$
l'approximation ainsi obtenue de $f'(x) \cdot v \in \mathbb{R}^{m}$ est d'ordre 2 si $g^{(3)}(0) \neq 0$ ou au moins d'ordre 4 sinon. À noter que ce schéma nécessite plus d'évaluations de la fonction $f$. On peut montrer comme précédemment que le meilleur $h$ est de l'ordre
$$
h_* \approx \sqrt[3]{\mathrm{eps}_\mathrm{mach}}.
$$
%% Cell type:markdown id: tags:
## Dérivées par différentiation complexe
Les formules des schémas avant et centrée sont sensibles aux calculs de la différence $\Delta f = f(x+h) - f(x)$ ou $\Delta f = f(x+h) - f(x-h)$. Pour remédier à ce problème, les [différences finies à pas complexe](https://dl.acm.org/doi/10.1145/838250.838251) ont été introduites.
Si on suppose que la fonction $g$ est holomorphe, c'est-à-dire dérivable au sens complexe,
on peut considérer un pas complexe $ih$. Un développement limité de $g$ en $0$ s'écrit
$$
f(x+ih v) = g(ih) = g(0) + ih g'(0) - \frac{h^2}{2} g^{(2)}(0) - i\frac{h^3}{6} g^{(3)}(0) + o(h^3),
$$
On considère alors l'approximation :
$$
f'(x) \cdot v = g'(0) \approx \frac{\mathrm{Im}(f(x+ihv))}{h}.
$$
On peut prouver que l'approximation ci-dessus est au moins d'ordre 2 et aussi démontrer que tout pas inférieur à $h_*$ est optimal, avec
$$
h_{*} \approx \sqrt{\mathrm{eps}_{\mathrm{mach}}}.
$$
**Remarque.** Utiliser en `Julia` la commande `imag` pour calculer la partie imaginaire d'un nombre complexe et la variable `im` pour représenter l'unité imaginaire.
%% Cell type:markdown id: tags:
## Dérivées par différentiation automatique via les nombres duaux
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) := f(a) + f'(a)\, b\, \varepsilon
$$
et
$$
g(a + b\, \varepsilon) := g(a) + g'(a)\, b\, \varepsilon.
$$
On a alors automatiquement 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$
Voici comment créer un nombre dual en `Julia` et récupérer les parties réelles et duales (avec ce que j'ai défini ci-dessous) :
```julia
using DualNumbers
# scalar case
d = 1 + 2ε # ou 1 + 2 * ε ou 1 + ε * 2
real(d) # 1
dual(d) # 2
# vector case
d = [1, 3] + [2, 4]ε # ou [1, 3] + [2, 4] * ε ou [1, 3] + ε * [2, 4] ou [1+2ε, 3+4ε]
real(d) # [1, 3]
dual(d) # [2, 4]
```
**Remarque.** On peut aussi utiliser le package `ForwardDiff` pour calculer des dérivées automatiquement. Il est plus performant que `DualNumbers`.
%% Cell type:markdown id: tags:
## Fonctions auxiliaires
%% Cell type:code id: tags:
``` julia
# available methods
methods = (:forward, :central, :complex, :dual, :forward_ad)
# type of x or its coordinates
function mytypeof(x::Union{T, Vector{<:T}}) where T<:AbstractFloat
return T
end
# default step value
function _step(x, v, method)
T = mytypeof(x)
eps_value = T isa AbstractFloat ? eps(T) : eps(1.)
if method == :forward
step = √(eps_value)
elseif method == :central
step = (eps_value)^(1/3)
elseif method == :complex
step = √(eps_value)
else
step = 0.0
end
step *= √(max(1., norm(x))) / √(max(1.0, norm(v)))
return step
end
# default method value
function _method()
return :forward
end;
# creation of dual number ε
import Base.*
*(e::Function, x::Union{Number, Vector{<:Number}}) = e(x)
*(x::Union{Number, Vector{<:Number}}, e::Function) = e(x)
ε(x=1) = begin
if x isa Number
return Dual.(0.0, x)
else
return Dual.(zeros(length(x)), x)
end
end
em = ε
dual(x::Union{Dual, Vector{<:Dual}}) = dualpart.(x)
real(x::Union{Dual, Vector{<:Dual}}) = realpart.(x);
```
%% Cell type:markdown id: tags:
## La méthode principale pour le calcul de dérivées
La fonction `derivative` ci-dessous calcule la dérivée directionnelle
$$
f'(x) \cdot v.
$$
%% Cell type:code id: tags:
``` julia
## TO COMPLETE
# compute directional derivative
function derivative(f, x, v; method=_method(), h=_step(x, v, method))
if method ∉ methods
error("Choose a valid method in ", methods)
end
if method == :forward
return f(x) # TO UPDATE
elseif method == :central
return f(x) # TO UPDATE
elseif method == :complex
return f(x) # TO UPDATE
elseif method == :dual
return f(x) # TO UPDATE
elseif method == :forward_ad
if x isa Number
return ForwardDiff.derivative(f, x)*v
else
return ForwardDiff.jacobian(f, x)*v
end
end
end;
```
%% Cell type:code id: tags:
``` julia
# function to print derivative values and errors
function print_derivatives(f, x, v, sol)
println("Hand derivative: ", sol, "\n")
for method ∈ methods
dfv = derivative(f, x, v, method=method)
println("Method: ", method)
println(" derivative: ", dfv)
@printf(" error: %e\n", norm(dfv - sol))
if method ∈ (:forward, :central, :complex)
step = _step(x, v, method)
println(" step: ", step)
@printf(" error/step: %e\n", norm(dfv - sol) / step)
end
println()
end
end;
```
%% Cell type:markdown id: tags:
## Exercice 1
1. Compléter la fonction `derivative` ci-dessus avec les méthodes de différences finies avant, centrée, par différentiation complexe et par différentiation automatique via les nombres duaux.
2. Exécuter le code ci-dessous et vérifier les résultats obtenus.
%% Cell type:code id: tags:
``` julia
# Scalar case
# check if the derivatives are correct
f(x) = cos(x)
x = π/4
v = 1.0
# solution
sol = -sin(x)*v
# print derivatives and errors for each method
print_derivatives(f, x, v, sol)
```
%% Cell type:code id: tags:
``` julia
# Vectorial case
# check if the derivatives are correct
f(x) = [0.5*(x[1]^2 + x[2]^2); x[1]*x[2]]
x = [1.0, 2.0]
v = [1.0, -1.0]
# solution
sol = [x[1]*v[1]+x[2]*v[2], x[1]*v[2]+x[2]*v[1]]
# print derivatives and errors for each method
print_derivatives(f, x, v, sol)
```
%% Cell type:markdown id: tags:
## Pas optimal
On se propose de tester pour la fonction $\cos$ aux points $x_0=\pi/3$, $x_1 = 10^6\times\pi/3$ et la fonction $\cos+10^{-8} \mathcal{N}(0, \, 1)$ au point $x_0=\pi/3$ l'erreur entre les différences finies et la dérivée au point considéré en fonction de $h$. On prendra $h=10^{-i}$ pour $i= \{1,\ldots,16\}$ et on tracera ces erreurs dans une échelle logarithmique (en `Julia`, avec le package `Plots` on utilise l'option `scale=:log10`).
## Exercice 2
- Visualiser les différentes erreurs en fonction de $h$ pour les différentes méthodes de calcul de dérivées. Commentaires.
- Modifier la précision de $x_0$ et $x_1$ en `Float32`. Commentaires.
%% Cell type:code id: tags:
``` julia
# affichage des erreurs en fonction de h
function plot_errors(steps, errors, h_star, title)
steps_save = steps
ymax = 10^10
# supprimer les erreurs nulles
non_nul_element = findall(!iszero, errors)
errors = errors[non_nul_element]
steps = steps[non_nul_element]
# Courbe des erreurs pour les differents steps en bleu
plt = plot((10.).^(-steps), errors, xscale=:log10, yscale=:log10, linecolor=:blue, lw=2, legend=false)
# régler xlims pour toujours avoir tous les steps de départ
plot!(plt, xlims=(10^(-maximum(steps_save)), 10^(-minimum(steps_save))))
# ylims toujours entre 10^-16 et ymax
plot!(plt, ylims=(10^(-16), ymax))
# Ligne verticale pour situer l'erreur optimale h* en rouge
plot!(plt,[h_star, h_star], [10^(-16), ymax], linecolor=:red, lw=1, linestyle=:dash)
# titre de la figure et xlabel
plot!(plt, xlabel = "h", title = title, legend=false, titlefontsize=10)
# ajouter des marges en bas de la figure pour mieux voir le xlabel
plot!(plt, bottom_margin = 5mm)
#
return plt
end;
```
%% Cell type:code id: tags:
``` julia
# Les differentes fonctions et la dérivée theorique
fun1(x) = cos(x)
fun2(x) = cos(x) + 1.e-8*randn()
dfun(x) = -sin(x);
```
%% Cell type:code id: tags:
``` julia
# method
method = :forward # TO PLAY WITH
h_star = √(eps(1.0)) # TO UPDATE ACCORDING TO THE METHOD
# Points pour lesquels on souhaite effectuer les tests
x0 = π/3
x1 = 1.e6*π/3
# steps pour faire les tests
steps = range(1, 16, 16)
# Initialisation des vecteurs d'erreur
err_x0 = zeros(length(steps))
err_x0p = zeros(length(steps))
err_x1 = zeros(length(steps))
# Calcul des erreurs
for i in 1:length(steps)
h = 10^(-steps[i])
err_x0[i] = abs(derivative(fun1, x0, 1.0, h=h, method=method) - (dfun(x0)))
err_x1[i] = abs(derivative(fun1, x1, 1.0, h=h, method=method) - (dfun(x1)))
err_x0p[i] = abs(derivative(fun2, x0, 1.0, h=h, method=method) - (dfun(x0)))
end
# Affichage des erreurs
p1 = plot_errors(steps, err_x0, h_star, "cos(x0)")
p2 = plot_errors(steps, err_x0p, h_star, "cos(x0) + perturbation")
p3 = plot_errors(steps, err_x1, h_star, "cos(x1)")
plot(p1, p2, p3, layout=(1,3), size=(850, 350))
```
%% Cell type:code id: tags:
``` julia
```
tp/direct-indirect.ipynb deleted 100644 → 0
+ 0
− 430
View file @ 578b17ca
%% Cell type:markdown id: tags:
[<img src="https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png" alt="N7" height="100"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)
# Méthode directe et tir simple structuré
- Date : 2024
- Durée approximative : 1.5 séance
Le but est de résoudre par du tir simple indirect, un problème de contrôle optimal dont le contrôle est discontinu. On se propose de résoudre dans un premier temps, le problème par une méthode directe, afin de déterminer la structure optimale et une bonne approximation de la solution pour faire converger la méthode de tir.
%% Cell type:code id: tags:
``` julia
using JuMP, Ipopt # pour la méthode directe
using DifferentialEquations, NLsolve, ForwardDiff # pour la méthode indirecte
using Plots # pour les graphiques
include("utils.jl"); # fonctions utilitaires
```
%% Cell type:markdown id: tags:
## Introduction
On considère le problème de contrôle optimal suivant :
$$
\left\{
\begin{array}{l}
\displaystyle \min_{x, u} \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = \displaystyle u(t), \quad |u(t)| \le 1, \quad t \in [0, t_f] \text{ p.p.}, \\[1.0em]
x(0) = 1, \quad x(t_f) = 1/2
\end{array}
\right.
$$
où $t_f$ = 2.
✏️ **Exercice 1.**
1. Appliquez le PMP. Que pouvez-vous dire du contrôle maximisant ?
2. Peut-on appliquer simplement la méthode de tir simple vu au TP précédent ?
%% Cell type:markdown id: tags:
## Méthode directe
Avant de définir la méthode directe, on propose une réécriture du problème. Notez que ce n'est pas une obligation en soi de la méthode mais cela simplifie les choses.
✏️ **Exercice 2.**
- Mettez le problème sous forme de Mayer (c.f. cours). Vous nommerez la nouvelle variable d'état associée au coût $c(\cdot)$.
%% Cell type:markdown id: tags:
**Description de la méthode directe.**
L'idée principale est de transformer le problème de contrôle optimal (de dimension infinie) en un problème d'optimisation de dimension finie.
Pour cela :
1. On définit une grille uniforme en temps $(t_1, \dots, t_{N+1})$ où $t_1 = 0$ et $t_{N+1} = t_f$ avec $\Delta t = (t_f - t_0)/N$ le pas de discrétisation.
2. On discrétise l'état, le contrôle et le coût sur cette grille et on note
$$
X = \{ (x_i, u_i, c_i) ~|~ i \in \{1, \dots, N+1\}\}
$$
l'ensemble des variables d'optimisation du problème discrétisé.
Si l'on note $(x^*(\cdot), u^*(\cdot), c^*(\cdot))$ la solution du problème de contrôle optimal et $\{ (x^*_i, u^*_i, c^*_i) ~|~ i \in \{1, \dots, N+1\}\}$ la solution du problème discrétisé, on s'attend à avoir
$$
x^*_i \approx x^*(t_i), \quad u^*_i \approx u^*(t_i), \quad c^*_i \approx c^*(t_i), \quad \forall i \in \{1, \dots, N+1\}.
$$
✏️ ️**Exercice 3.** Définissez pour le problème discrétisé :
1. l'objectif à minimiser en fonction d'une ou plusieurs composantes de $X$ bien choisies.
2. les contraintes d'inégalités associées au contrôle.
3. les contraintes initiales et finales associées à la variable d'état.
4. les contraintes de dynamique sur l'état et le coût, en utilisant le schéma d'intégration de Crank-Nicolson (ou [règle des Trapèzes](https://fr.wikipedia.org/wiki/M%C3%A9thode_des_trap%C3%A8zes)).
%% Cell type:markdown id: tags:
✏️ ️**Exercice 4.**
- Modifiez le code suivant afin de résoudre le problème par la méthode directe que vous venez de décrire.
**Remarque.** Vous pouvez vous inspirer de cet [exemple](https://ct.gitlabpages.inria.fr/gallery/goddard-j/goddard.html) pour le code. Notez que dans cet exemple, la fonction `@NLexpressions` est utilisée mais n'est pas nécessaire ici donc vous pouvez ou non l'utiliser.
%% Cell type:code id: tags:
``` julia
# Création du modèle JuMP, Utilisation de Ipopt comme solver
sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5))
set_optimizer_attribute(sys,"tol",1e-8)
set_optimizer_attribute(sys,"constr_viol_tol",1e-6)
set_optimizer_attribute(sys,"max_iter",1000)
#######
####### DEBUT A MODIFIER
#######
####### Les ... sont à remplacer !
#######
####### Attention, on écrit le problème sous la forme d'un problème de Mayer
#######
# Paramètres
t0 = 0 # temps initial
tf = 0 # temps final
c0 = 0 # coût initial
x0 = 0 # état initial
xf = 0 # état final
N = 50 # taille de la grille
Δt = (tf-t0)/N # pas de temps
@variables(sys, begin
... # coût
... # état
... # control
end)
# Objectif
@objective(sys, Min, ...)
# Conditions initiales et finales
@constraints(sys, begin
con_c0, ... # contraint sur le coût initial
con_x0, ... # contraint sur l'état initial
con_xf, ... # contraint sur l'état final
end)
# Contraintes de dynamique avec le schéma de Crank-Nicolson
@NLconstraints(sys, begin
con_dc[j=1:N], ... # contraint différentielle sur le coût
con_dx[j=1:N], ... # contraint différentielle sur l'état
end);
#######
####### FIN A MODIFIER
#######
```
%% Cell type:code id: tags:
``` julia
# Solve for the control and state
println("Solving...")
optimize!(sys)
println()
# Display results
if termination_status(sys) == MOI.OPTIMAL
println(" Solution is optimal")
elseif termination_status(sys) == MOI.LOCALLY_SOLVED
println(" (Local) solution found")
elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys)
println(" Solution is suboptimal due to a time limit, but a primal solution is available")
else
error(" The model was not solved correctly.")
end
println(" objective value = ", objective_value(sys))
println()
# Retrieves values (including duals)
c = value.(c)[:]
x = value.(x)[:]
u = value.(u)[:]
t = (0:N) * value.(Δt)
pc0 = dual(con_c0)
px0 = dual(con_x0)
pxf = dual(con_xf)
if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end
if(px0*dual(con_dx[1])<0); px0 = -px0; end
if(pxf*dual(con_dx[N])<0); pxf = -pxf; end
if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle
sign = -1.0
else
sign = 1.0
end
pc = [ dual(con_dc[i]) for i in 1:N ]
px = [ dual(con_dx[i]) for i in 1:N ]
pc = sign * [pc0; pc[1:N]]
px = sign * [px0; (px[1:N-1]+px[2:N])/2; pxf]; # We add the multiplier from the limit conditions
```
%% Cell type:code id: tags:
``` julia
x_plot = plot(t, x, xlabel = "t", ylabel = "x", legend = false)
u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre)
px_plot = plot(t, px, xlabel = "t", ylabel = "p", legend = false)
display(plot(x_plot, px_plot, layout = (1,2), size=(800,300)))
display(u_plot)
```
%% Cell type:markdown id: tags:
✏️ ️**Exercice 5.**
1. Commentez les résultats.
2. Modifiez la tolérance de l'optimiseur ainsi que le nombre de points de discrétisation. Commentaires.
3. Décrivez la structure du contrôle optimal en fonction du temps.
%% Cell type:markdown id: tags:
## Méthode indirecte
D'après la condition de maximisation du hamiltonien et la structure optimale de la solution, nous avons besoin de définir 3 fonctions permettant de calculer les flots des systèmes hamiltoniens associés aux contrôles +1, -1 et au contrôle dit singulier qui apparaît lorsque la fonction de commutation reste nulle sur un intervalle de temps non réduit à un singleton.
✏️ ️️**Exercice 6.**
1. Donner la fonction de commutation.
2. En dérivant deux fois par rapport au temps la fonction de commutations, trouver une condition vérifiée par l'état et une condition vérifiée par le contrôle le long d'un arc singulier, c-à-d le long d'un arc où la fonction de commutation est constante égale à 0.
3. Remplir le code ci-dessous: il faut fournir les 3 contrôles pour définir les flots, puis écrire la fonction de tir.
**Remarque.** La fonction de tir à 3 variables inconnues. Il faut donc trouver 3 équations. La condition terminale en est une. L'annulation de la fonction de commutation au premier temps de commutation en est une autre. La question 2 aide à trouver la dernière condition.
%% Cell type:code id: tags:
``` julia
# système pseudo-hamiltonien
function hv(x, p, u)
return [u, 2x]
end
# systèmes hamiltoniens
hv_min(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
hv_max(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
hv_sin(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
# flots
fmin = Flow(hv_min)
fmax = Flow(hv_max)
fsin = Flow(hv_sin)
# fonction de tir
function shoot(p0, t1, t2)
# integration
x1, p1 = fmin(t0, x0, p0, t1) # x1, p1 sont des scalaires
x2, p2 = fsin(t1, x1, p1, t2)
x3, p3 = fmax(t2, x2, p2, tf)
# conditions
s = zeros(eltype(p0), 3)
s[1] = NaN ## REMPLACER NaN par la bonne expression
s[2] = NaN ## REMPLACER NaN par la bonne expression
s[3] = NaN ## REMPLACER NaN par la bonne expression
return s
end;
```
%% Cell type:markdown id: tags:
4. Trouver (remplir le code ci-dessous) à partir de la solution de la méthode directe, une bonne initialisation pour la fonction de tir.
%% Cell type:code id: tags:
``` julia
# itéré initiale pour la méthode indirecte de tir
p0 = NaN ## REMPLACER NaN par une bonne valeur
t1 = NaN ## REMPLACER NaN par une bonne valeur
t2 = NaN ## REMPLACER NaN par une bonne valeur
#
y = [ p0 ; t1 ; t2]
println("Itéré initial:\n", y)
```
%% Output
Itéré initial:
[NaN, NaN, NaN]
%% Cell type:code id: tags:
``` julia
# fonction de tir et sa jacobienne
foo(y) = shoot(y...)
jfoo(y) = ForwardDiff.jacobian(foo, y)
# Résolution de shoot(p0, t1, t2) = 0.
nl_sol = nlsolve(foo, jfoo, y; xtol=1e-8, method=:trust_region, show_trace=true);
# Retrieves solution
if converged(nl_sol)
p0 = nl_sol.zero[1]
t1 = nl_sol.zero[2]
t2 = nl_sol.zero[3];
println("\nFinal solution:\n", nl_sol.zero)
else
error("Not converged")
end
```
%% Cell type:code id: tags:
``` julia
# affichage de la solution
ode_sol = fmin((t0, t1), x0, p0)
tt0 = ode_sol.t
xx0 = [ ode_sol[1, j] for j in 1:size(tt0, 1) ]
pp0 = [ ode_sol[2, j] for j in 1:size(tt0, 1) ]
uu0 = -ones(size(tt0, 1))
ode_sol = fsin((t1, t2), xx0[end], pp0[end])
tt1 = ode_sol.t
xx1 = [ ode_sol[1, j] for j in 1:size(tt1, 1) ]
pp1 = [ ode_sol[2, j] for j in 1:size(tt1, 1) ]
uu1 = zeros(size(tt1, 1))
ode_sol = fmax((t2, tf), xx1[end], pp1[end])
tt2 = ode_sol.t
xx2 = [ ode_sol[1, j] for j in 1:size(tt2, 1) ]
pp2 = [ ode_sol[2, j] for j in 1:size(tt2, 1) ]
uu2 = ones(size(tt2, 1))
t_shoot = [ tt0 ; tt1 ; tt2 ]
x_shoot = [ xx0 ; xx1 ; xx2 ]
p_shoot = [ pp0 ; pp1 ; pp2 ]
u_shoot = [ uu0 ; uu1 ; uu2 ]
x_plot = plot(t_shoot, x_shoot, xlabel = "t", ylabel = "x", legend = false)
u_plot = plot(t_shoot, u_shoot, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre)
px_plot = plot(t_shoot, p_shoot, xlabel = "t", ylabel = "px", legend = false)
display(plot(x_plot, px_plot, layout = (1,2), size=(800, 300)))
display(u_plot)
```
%% Cell type:markdown id: tags:
5. Comparer avec l'affichage de la solution de la méthode directe.
tp/rk-implicites.ipynb deleted 100644 → 0
+ 0
− 467
View file @ 578b17ca
%% Cell type:markdown id: tags:
<img src="https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png" alt="N7" height="80"/>
<img src="https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png" alt="INSA" height="80"/>
# Méthodes de Runge-Kutta implicites - TP Projet
- Date : 2023-2024
- Durée approximative : inconnue
%% Cell type:markdown id: tags:
**Nom** : "Entrez votre nom ici"
**Prénom** : "Entrez votre prénom ici"
%% Cell type:markdown id: tags:
## Rendu et consignes
Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.
- **Date limite du rendu : mercredi 15/11/2023 à 23h59.** Attention, à partir de 24h, 2 points est enlevé de la note finale toutes les 30 minutes.
- **Attention :** Le fichier doit être renommé de la façon suivante : `rk_implicites_NOM_Prenom.ipynb`. 4 points enlevés si le nom du fichier n'est pas respecté.
- **Documents autorisés :** vous pouvez utiliser votre polycopié et les TPs précédents.
%% Cell type:code id: tags:
``` julia
a = (1, 2)
b = (3, 4)
a.-b
```
%% Output
%% Cell type:code id: tags:
``` julia
f(x, y) = x+y
```
%% Output
%% Cell type:code id: tags:
``` julia
f(a...)
```
%% Output
%% Cell type:code id: tags:
``` julia
v = [a...]
```
%% Output
%% Cell type:code id: tags:
``` julia
using LinearAlgebra
k1 = [1, 2]
k2 = [3, 4]
ks = (k1, k2)
G(k1, k2) = (k1, k2)
norm(ks.-G(ks...))
```
%% Output
%% Cell type:markdown id: tags:
## Introduction
Nous allons dans ce TP, implémenter quelques méthodes de Runge-Kutta **implicites** (voir **polycopié Section 8.2**) et étudier leur convergence. On considère un pas de temps $h$ uniforme. Une méthode à un pas implicite est convergente si pour toute solution $x(\cdot, t_0, x_0)$ d'un problème de Cauchy, la suite approximante ${(x_i)}_i$ donnée par la méthode à un pas implicite vérifie
$$
\max_{1 \le i \le N}\, \|{x(t_i, t_0, x_0) - x_i}\| \to 0
\quad\text{quand}\quad h \to 0.
$$
Si la convergence est d'ordre $p$ alors il existe une constante $C \ge 0$ telle que l'**erreur globale** $E$ vérifie
$$
E := \max_{1 \le i \le N}\, \|{x(t_i, t_0, x_0) - x_i}\| \le C\, h^p.
$$
Faisons l'hypothèse que $E = M\, h^p$ pour un certain $M \ge 0$. En passant au logarithme, on obtient
$$
\log(E) = \log(M) + p\, \log(h).
$$
Nous en déduisons que si on trace $\log(E)$ en fonction de $\log(h)$, on doit obtenir une droite de pente $p$. C'est ce que nous allons vérifier dans ce TP.
%% Cell type:code id: tags:
``` julia
# activate local project
using Pkg
Pkg.activate(".")
# load packages
using LinearAlgebra
using Plots
using Plots.PlotMeasures
using Polynomials
#
px = PlotMeasures.px;
```
%% Cell type:code id: tags:
``` julia
# Fonctons auxiliaires
# VOUS POUVEZ METTRE VOS FONCTIONS AUXILIAIRES ICI
```
%% Cell type:markdown id: tags:
## L'exemple d'étude
On s'intéresse (pour les exercices 1, 2 et 3) au problème de Cauchy
$$
\dot{x}(t) = (1-2t) x(t), \quad x(0) = x_0 = 1
$$
sur l'intervalle $[0, 3]$.
%% Cell type:code id: tags:
``` julia
# Définition du problème de Cauchy
f(t, x) = (1-2t)x # Second membre f(t, x)
x0 = 1.0 # Condition initiale
tspan = (0.0, 3.0); # Intervalle de temps
```
%% Cell type:code id: tags:
``` julia
# Solution analytique
function sol(t)
return exp(t-t^2)
end;
```
%% Cell type:code id: tags:
``` julia
# Estimation de la constante de Lipschitz de f sur [0, 3]
# Voir Théorème 8.2.2 pour l'utilité de cette estimation
function dfx(t)
return 1-2t
end
L = maximum(abs.(dfx.(range(0, stop=3, length=1000))))
```
%% Cell type:markdown id: tags:
## La méthode d'Euler implicite
La méthode d'Euler implicite est donnée par :
$$
\left\{
\begin{array}{l}
x_{n+1} = x_n + h f(t_n + h, x_{n+1}) = G(x_{n+1}) \\
x_0 = x(t_0)
\end{array}
\right.
$$
### Exercice 1
1. Implémenter la méthode d'Euler implicite avec le point fixe (penser à voir le polycopié Section 8.2).
2. Pourquoi si $h \ge 0.2$, l'algorithme du point fixe peut s'arrêter au bout du nombre d'itérations max et donc ne pas converger pour la méthode d'Euler implicite ?
3. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.
4. Tracer l'erreur globale de la méthode d'Euler implicite en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$.
**Attention** : pour l'algorithme du point fixe, faites attention aux critères d'arrêts (il y en a 2) ! Voir votre polycopié Section 8.2. Vous fixerez la valeur de la tolérance à $10^{-6}$ et le nombre maximum d'itérations à $1000$.
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
# F(y) = y - G(y)
G(y) = (y + sin(y))/3
y = 1
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
y = G(y)
```
%% Output
%% Cell type:code id: tags:
``` julia
F(y) = y - G(y)
F(y)
```
%% Output
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## La méthode des trapèzes
La méthode des trapèzes est donnée par le tableau de Butcher :
$$
\begin{array}{c | c c}
0 & 0 & 0 \\[0.2em]
1 & 1/2 & 1/2 \\[0.2em]
\hline
& 1/2 & 1/2 \\
\end{array}
$$
### Exercice 2
1. Implémenter la méthode des trapèzes avec le point fixe.
2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.
3. Tracer l'erreur globale de la méthode des trapèzes. Quel est l'ordre de convergence de la méthode des trapèzes ?
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
## La méthode de Gauss à 2 étages
La méthode de Gauss à 2 étages est donnée par le tableau de Butcher :
$$
\begin{array}{c | c c}
1/2 - \sqrt{3}/6 & 1/4 & 1/4 - \sqrt{3}/6 \\[0.2em]
1/2 + \sqrt{3}/6 & 1/4 + \sqrt{3}/6 & 1/4 \\[0.2em]
\hline
& 1/2 & 1/2 \\
\end{array}
$$
### Exercice 3
1. Implémenter la méthode de Gauss à 2 étages avec le point fixe.
2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.
3. Tracer l'erreur globale de la méthode de Gauss à 2 étages. Quel est l'ordre de convergence de la méthode de Gauss à 2 étages ?
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
%% Cell type:markdown id: tags:
## Un autre exemple
On considère à partir de maintenant l'équation différentielle en dimension 2 :
$$
\dot{x}_1(t) = x_2(t), \quad \dot{x}_2(t) = - x_1(t).
$$
On peut montrer facilement que la norme de $x(t) = (x_1(t), x_2(t))$ est constante le long des solutions :
$$
\frac{\mathrm{d}}{\mathrm{d} t} \|x(t)\|^2 = 2\, \left( x(t) \,|\, \dot{x}(t) \right) = 2 \left( x_1(t) x_2(t) - x_2(t) x_1(t) \right) = 0.
$$
### Exercice 4
On considère le problème de Cauchy associé de condition initiale $x_0 = (1, 0)$.
1. Afficher l'approximation de la solution sur $[0, 10]$ pour les méthodes :
- Euler explicite ;
- Euler implicite ;
- Trapèzes ;
- Gauss à 2 étages.
2. Commentaires.
**Attention :** vous ferez un affichage dans le plan $(x_1, x_2)$. Vous fixerez le nombre de pas à $N=100$.
%% Cell type:code id: tags:
``` julia
# AJOUTER VOTRE CODE ICI
```
tp/runge-kutta.ipynb deleted 100644 → 0
+ 0
− 771
View file @ 578b17ca
This diff is collapsed.
tp/simple-shooting.ipynb deleted 100644 → 0
+ 0
− 368
View file @ 578b17ca
%% Cell type:markdown id:efcf7482 tags:
[<img src="https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png" alt="N7" height="100"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)
# Tir simple indirect
- Date : 2023
- Durée approximative : 3/4 de séance
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.
%% Cell type:code id:1e3057d1 tags:
``` julia
# Packages
using DifferentialEquations
using Plots
using NLsolve
using ForwardDiff
using LinearAlgebra
include("utils.jl"); # plot_traj!, plot_flow!, Flow
```
%% 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.
**Remarque.** Utiliser le TP précédent pour voir comment répondre à la question. Utiliser notamment les fonctions `plot_traj!` et `plot_flot!` du fichier `utils.jl`.
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. Le point $z(t_f, x_0, p_0)$ est la solution du système hamiltonien au temps $t_f$ qui part au temps initial $t_0=0$ du point $(x_0, p_0)$. Il est calculé en `Julia` à l'aide de la fonction $f$ ci-dessous, obtenue par `f = Flow(Hv)`.
3. Afficher la fonction de tir pour $p_0 \in [0, 1]$.
4. Résoudre l'équation de tir $S(p_0) = 0$ et tracer l'extrémale dans le plan de phase, cf. question 1.
**Remarque.** Donner une bonne initialisation au solveur de Newton pour la résolution de l'équation de tir.
%% 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; # correspond à xf dans la description mathématique du problème
```
%% Cell type:code id:497d73ec-5b23-459e-8c32-ca13add12f5f tags:
``` julia
## Question 1:
# ------------------------------------
# conseils :
# - revoir la fin du TP précédent, ie le TP ode
# - réutiliser plot_traj! et plot_flow! comme au TP ode, cf utils.jl
#
# ------------------------------------
# système hamiltonien sous la forme F(z), z = (x, p)
#
# à 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: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 = ẋ, dp = ṗ
"""
dx = 0
dp = 0
return dx, dp
end
# flot du système hamiltonien: voir utils.jl
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 finaliser l'affichage
```
%% 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: contrainte du solveur
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, cf. f = Flow(Hv). Voir aussi utils.jl.
#
plt # A REMPLACER: attention, il faut remettre les xlims et ylims pour un affichage propre
```
%% 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
f = 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 = f((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))
```
tp/space.jl deleted 100644 → 0
+ 0
− 849
View file @ 578b17ca
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment