diff --git a/README.md b/README.md
index ed3594932c3c330a7e542d92aac8bcc8936b55d7..ecb43417209c421d9ddd6cdae5e74a811865b9e7 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,9 @@ Cours de contrôle optimal de l'ENSEEIHT pour le département Sciences du Numér
**TP**
-* [TP1 : tp/ode.ipynb](tp/ode.ipynb) - Intégration numérique et systèmes dynamiques
+* [Sujet 1 : tp/ode.ipynb](tp/ode.ipynb) - Intégration numérique et systèmes dynamiques
+* [Sujet 2 : tp/simple-shooting.ipynb](tp/simple-shooting.ipynb) - Tir simple indirect
+* [Sujet 3 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul numérique de dérivées
**Notes de cours supplémentaires - ressources**
diff --git a/tp/derivative.ipynb b/tp/derivative.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..f4869bb4a8c96caf4fc5c792afa7fa5da3d52cf2
--- /dev/null
+++ b/tp/derivative.ipynb
@@ -0,0 +1,495 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "[<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)\n",
+ "\n",
+ "# Calcul de dérivées\n",
+ "\n",
+ "- Date : 2023-2024\n",
+ "- Durée approximative : 1h15\n",
+ "\n",
+ "## Introduction\n",
+ "\n",
+ "Il existe plusieurs façon de calculer une dérivée sur un calculateur : \n",
+ "\n",
+ "- par différences finies;\n",
+ "- par différentiation complexe;\n",
+ "- en utilisant la différentiation automatique;\n",
+ "- en utilisant le calcul formel et un générateur de code.\n",
+ "\n",
+ "Nous allons étudier ici quelques cas.\n",
+ "\n",
+ "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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# load packages\n",
+ "using DualNumbers\n",
+ "using DifferentialEquations\n",
+ "using ForwardDiff\n",
+ "using LinearAlgebra\n",
+ "using Plots\n",
+ "using Plots.Measures\n",
+ "using Printf"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Dérivées par différences finies avant\n",
+ "\n",
+ "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 :\n",
+ "$$\n",
+ " g(h) = \\sum_{i=0}^{n} \\frac{h^i}{i!} g^{(i)}(0) + R_n(h), \\quad R_n(h) = o(h^n),\n",
+ "$$\n",
+ "ou d'après Taylor-Lagrange, \n",
+ "$$\n",
+ " \\| R_n(h) \\| \\leq \\frac{M_n h^{n+1}}{(n+1)!},\n",
+ "$$\n",
+ "de même, \n",
+ "$$\n",
+ " g(-h) = \\sum_{i=0}^{n} \\frac{(-h)^i}{i!} g^{(i)}(0) + R_n(h).\n",
+ "$$\n",
+ "\n",
+ "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 : \n",
+ "$$\n",
+ " \\frac{f(x+hv) - f(x)}{h} = \n",
+ " \\frac{g(h)-g(0)}{h} = g'(0) + \\frac{h}{2} g^{2}(0) + \\frac{h^2}{6} g^{(3)}(0) + o(h^2).\n",
+ "$$\n",
+ "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. \n",
+ "\n",
+ "**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.\n",
+ "\n",
+ "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 : \n",
+ "$$\n",
+ " \\left\\| \\mathrm{num}(g,h) - g(h) \\right\\| := \\| e_h\\| \\leq \\mathrm{eps}_\\mathrm{mach} L_f,\n",
+ "$$\n",
+ "ou $L_f$ est une constante qui dépend de la valeur de $f$ sur le domaine d'intérêt. Ainsi on a : \n",
+ "\\begin{align*}\n",
+ " \\left\\| \\frac{\\mathrm{num}(g,h) - \\mathrm{num}(g,0)}{h} - g'(0) \\right\\|\n",
+ " &= \\left\\| \\frac{g(h) + e_h - g(0) - e_0}{h} - g'(0) \\right\\|, \\\\[1em]\n",
+ " &= \\left\\| \\frac{R_1(h)}{h} + \\frac{e_h - e_0}{h} \\right\\|, \\\\[1em]\n",
+ " &\\leq \\left\\| \\frac{R_1(h)}{ h} \\right\\| + \\left\\| \\frac{e_h - e_0}{h} \\right\\|, \\\\[1em]\n",
+ " & \\leq \n",
+ " \\underbrace{ \\frac{M_1 h}{2}}_{{\\text{Erreur d'approximation}}} + \n",
+ " \\underbrace{2 \\frac{\\mathrm{eps}_\\mathrm{mach}L_f}{h}}_{{\\text{Erreur numérique}}}.\n",
+ "\\end{align*} \n",
+ "Le majorant trouvé atteint son minimum en \n",
+ "$$\n",
+ " h_{*} = 2 \\sqrt{\\frac{\\mathrm{eps}_\\mathrm{mach} L_f }{M_1}}.\n",
+ "$$\n",
+ "\n",
+ "En considérant que $L_f \\simeq M_1$, alors le choix se révélant le plus optimal est \n",
+ "$$\n",
+ " h_{*} \\approx \\sqrt{\\mathrm{eps}_\\mathrm{mach}}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Dérivées par différences finies centrées\n",
+ "\n",
+ "On peut utiliser un schéma de différences finies centrée pour calculer la dérviée de $g$. \n",
+ "$$\n",
+ " \\frac{f(x+hv) - f(x-hv)}{2h} = \\frac{g(h) - g(-h)}{2h} = \n",
+ " g'(0) + g^{(3)}(0) \\frac{h^2}{6} + \\mathcal{O}(h^4),\n",
+ "$$\n",
+ "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 \n",
+ "$$\n",
+ " h_* \\approx \\sqrt[3]{\\mathrm{eps}_\\mathrm{mach}}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Dérivées par différentiation complexe\n",
+ "\n",
+ "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.\n",
+ "Si on suppose que la fonction $g$ est holomorphe, c'est-à-dire dérivable au sens complexe,\n",
+ "on peut considérer un pas complexe $ih$. Un développement limité de $g$ en $0$ s'écrit\n",
+ "\n",
+ "$$\n",
+ " 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),\n",
+ "$$\n",
+ "\n",
+ "On considère alors l'approximation : \n",
+ "\n",
+ "$$\n",
+ " f'(x) \\cdot v = g'(0) \\approx \\frac{\\mathrm{Im}(f(x+ihv))}{h}.\n",
+ "$$\n",
+ "\n",
+ "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 \n",
+ "$$\n",
+ " h_{*} \\approx \\sqrt{\\mathrm{eps}_{\\mathrm{mach}}}.\n",
+ "$$\n",
+ "\n",
+ "**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",
+ "metadata": {},
+ "source": [
+ "## Dérivées par différentiation automatique via les nombres duaux\n",
+ "\n",
+ "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.\n",
+ "\n",
+ "Soit deux fonctions $f, g \\colon \\mathbb{R} \\to \\mathbb{R}$ dérivables, de dérivées respectives $f'$ et $g'$. On pose\n",
+ "\n",
+ "$$\n",
+ "f(a + b\\, \\varepsilon) := f(a) + f'(a)\\, b\\, \\varepsilon\n",
+ "$$\n",
+ "\n",
+ "et\n",
+ "\n",
+ "$$\n",
+ "g(a + b\\, \\varepsilon) := g(a) + g'(a)\\, b\\, \\varepsilon.\n",
+ "$$\n",
+ "\n",
+ "On a alors automatiquement les propriétés suivantes. Posons $d = x + \\varepsilon$, alors :\n",
+ "\n",
+ "- $(f + g)(d) = (f+g)(x) + (f+g)'(x) \\, \\varepsilon$\n",
+ "- $(fg)(d) = (fg)(x) + (fg)'(x) \\, \\varepsilon$\n",
+ "- $(g \\circ f)(d) = (g \\circ f)(x) + (g \\circ f)'(x) \\, \\varepsilon$\n",
+ "\n",
+ "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) :\n",
+ "\n",
+ "```julia\n",
+ "using DualNumbers\n",
+ "\n",
+ "# scalar case\n",
+ "d = 1 + 2ε # ou 1 + 2 * ε ou 1 + ε * 2\n",
+ "real(d) # 1\n",
+ "dual(d) # 2\n",
+ "\n",
+ "# vector case\n",
+ "d = [1, 3] + [2, 4]ε # ou [1, 3] + [2, 4] * ε ou [1, 3] + ε * [2, 4] ou [1+2ε, 3+4ε]\n",
+ "real(d) # [1, 3]\n",
+ "dual(d) # [2, 4]\n",
+ "```\n",
+ "\n",
+ "**Remarque.** On peut aussi utiliser le package `ForwardDiff` pour calculer des dérivées automatiquement. Il est plus performant que `DualNumbers`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Fonctions auxiliaires"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# available methods\n",
+ "methods = (:forward, :central, :complex, :dual, :forward_ad)\n",
+ "\n",
+ "# type of x or its coordinates\n",
+ "function mytypeof(x::Union{T, Vector{<:T}}) where T<:AbstractFloat\n",
+ " return T\n",
+ "end\n",
+ "\n",
+ "# default step value\n",
+ "function _step(x, v, method)\n",
+ " T = mytypeof(x)\n",
+ " eps_value = T isa AbstractFloat ? eps(T) : eps(1.)\n",
+ " if method == :forward\n",
+ " step = √(eps_value)\n",
+ " elseif method == :central\n",
+ " step = (eps_value)^(1/3)\n",
+ " elseif method == :complex\n",
+ " step = √(eps_value)\n",
+ " else\n",
+ " step = 0.0\n",
+ " end\n",
+ " step *= √(max(1., norm(x))) / √(max(1.0, norm(v)))\n",
+ " return step\n",
+ "end\n",
+ "\n",
+ "# default method value\n",
+ "function _method()\n",
+ " return :forward \n",
+ "end;\n",
+ "\n",
+ "# creation of dual number ε\n",
+ "import Base.*\n",
+ "*(e::Function, x::Union{Number, Vector{<:Number}}) = e(x)\n",
+ "*(x::Union{Number, Vector{<:Number}}, e::Function) = e(x)\n",
+ "ε(x=1) = begin \n",
+ " if x isa Number\n",
+ " return Dual.(0.0, x)\n",
+ " else\n",
+ " return Dual.(zeros(length(x)), x)\n",
+ " end\n",
+ "end\n",
+ "em = ε\n",
+ "dual(x::Union{Dual, Vector{<:Dual}}) = dualpart.(x)\n",
+ "real(x::Union{Dual, Vector{<:Dual}}) = realpart.(x);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode principale pour le calcul de dérivées\n",
+ "\n",
+ "La fonction `derivative` ci-dessous calcule la dérivée directionnelle\n",
+ "\n",
+ "$$\n",
+ " f'(x) \\cdot v.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## TO COMPLETE\n",
+ "\n",
+ "# compute directional derivative\n",
+ "function derivative(f, x, v; method=_method(), h=_step(x, v, method))\n",
+ " if method ∉ methods \n",
+ " error(\"Choose a valid method in \", methods)\n",
+ " end\n",
+ " if method == :forward\n",
+ " return f(x) # TO UPDATE\n",
+ " elseif method == :central\n",
+ " return f(x) # TO UPDATE\n",
+ " elseif method == :complex\n",
+ " return f(x) # TO UPDATE\n",
+ " elseif method == :dual \n",
+ " return f(x) # TO UPDATE\n",
+ " elseif method == :forward_ad\n",
+ " if x isa Number\n",
+ " return ForwardDiff.derivative(f, x)*v\n",
+ " else\n",
+ " return ForwardDiff.jacobian(f, x)*v\n",
+ " end\n",
+ " end\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# function to print derivative values and errors\n",
+ "function print_derivatives(f, x, v, sol)\n",
+ "\n",
+ " println(\"Hand derivative: \", sol, \"\\n\")\n",
+ "\n",
+ " for method ∈ methods\n",
+ " dfv = derivative(f, x, v, method=method)\n",
+ " println(\"Method: \", method)\n",
+ " println(\" derivative: \", dfv)\n",
+ " @printf(\" error: %e\\n\", norm(dfv - sol))\n",
+ " if method ∈ (:forward, :central, :complex)\n",
+ " step = _step(x, v, method)\n",
+ " println(\" step: \", step)\n",
+ " @printf(\" error/step: %e\\n\", norm(dfv - sol) / step)\n",
+ " end\n",
+ " println()\n",
+ " end\n",
+ "\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercice 1\n",
+ "\n",
+ "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.\n",
+ "2. Exécuter le code ci-dessous et vérifier les résultats obtenus."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Scalar case\n",
+ "\n",
+ "# check if the derivatives are correct\n",
+ "f(x) = cos(x)\n",
+ "x = π/4\n",
+ "v = 1.0\n",
+ "\n",
+ "# solution\n",
+ "sol = -sin(x)*v\n",
+ "\n",
+ "# print derivatives and errors for each method\n",
+ "print_derivatives(f, x, v, sol)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Vectorial case\n",
+ "\n",
+ "# check if the derivatives are correct\n",
+ "f(x) = [0.5*(x[1]^2 + x[2]^2); x[1]*x[2]]\n",
+ "x = [1.0, 2.0]\n",
+ "v = [1.0, -1.0]\n",
+ "\n",
+ "# solution\n",
+ "sol = [x[1]*v[1]+x[2]*v[2], x[1]*v[2]+x[2]*v[1]]\n",
+ "\n",
+ "# print derivatives and errors for each method\n",
+ "print_derivatives(f, x, v, sol)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Pas optimal\n",
+ "\n",
+ "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`).\n",
+ "\n",
+ "## Exercice 2\n",
+ "\n",
+ "- Visualiser les différentes erreurs en fonction de $h$ pour les différentes méthodes de calcul de dérivées. Commentaires.\n",
+ "- Modifier la précision de $x_0$ et $x_1$ en `Float32`. Commentaires."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# affichage des erreurs en fonction de h\n",
+ "function plot_errors(steps, errors, h_star, title)\n",
+ "\n",
+ " steps_save = steps\n",
+ " ymax = 10^10\n",
+ "\n",
+ " # supprimer les erreurs nulles\n",
+ " non_nul_element = findall(!iszero, errors) \n",
+ " errors = errors[non_nul_element]\n",
+ " steps = steps[non_nul_element]\n",
+ "\n",
+ " # Courbe des erreurs pour les differents steps en bleu\n",
+ " plt = plot((10.).^(-steps), errors, xscale=:log10, yscale=:log10, linecolor=:blue, lw=2, legend=false)\n",
+ "\n",
+ " # régler xlims pour toujours avoir tous les steps de départ\n",
+ " plot!(plt, xlims=(10^(-maximum(steps_save)), 10^(-minimum(steps_save))))\n",
+ "\n",
+ " # ylims toujours entre 10^-16 et ymax\n",
+ " plot!(plt, ylims=(10^(-16), ymax))\n",
+ "\n",
+ " # Ligne verticale pour situer l'erreur optimale h* en rouge\n",
+ " plot!(plt,[h_star, h_star], [10^(-16), ymax], linecolor=:red, lw=1, linestyle=:dash)\n",
+ "\n",
+ " # titre de la figure et xlabel\n",
+ " plot!(plt, xlabel = \"h\", title = title, legend=false, titlefontsize=10)\n",
+ "\n",
+ " # ajouter des marges en bas de la figure pour mieux voir le xlabel \n",
+ " plot!(plt, bottom_margin = 5mm)\n",
+ "\n",
+ " #\n",
+ " return plt\n",
+ "\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Les differentes fonctions et la dérivée theorique\n",
+ "fun1(x) = cos(x)\n",
+ "fun2(x) = cos(x) + 1.e-8*randn()\n",
+ "dfun(x) = -sin(x);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# method \n",
+ "method = :forward # TO PLAY WITH\n",
+ "h_star = √(eps(1.0)) # TO UPDATE ACCORDING TO THE METHOD\n",
+ "\n",
+ "# Points pour lesquels on souhaite effectuer les tests\n",
+ "x0 = π/3\n",
+ "x1 = 1.e6*π/3\n",
+ "\n",
+ "# steps pour faire les tests\n",
+ "steps = range(1, 16, 16)\n",
+ "\n",
+ "# Initialisation des vecteurs d'erreur\n",
+ "err_x0 = zeros(length(steps))\n",
+ "err_x0p = zeros(length(steps))\n",
+ "err_x1 = zeros(length(steps))\n",
+ "\n",
+ "# Calcul des erreurs\n",
+ "for i in 1:length(steps)\n",
+ " h = 10^(-steps[i])\n",
+ " err_x0[i] = abs(derivative(fun1, x0, 1.0, h=h, method=method) - (dfun(x0)))\n",
+ " err_x1[i] = abs(derivative(fun1, x1, 1.0, h=h, method=method) - (dfun(x1)))\n",
+ " err_x0p[i] = abs(derivative(fun2, x0, 1.0, h=h, method=method) - (dfun(x0)))\n",
+ "end\n",
+ "\n",
+ "# Affichage des erreurs\n",
+ "p1 = plot_errors(steps, err_x0, h_star, \"cos(x0)\")\n",
+ "p2 = plot_errors(steps, err_x0p, h_star, \"cos(x0) + perturbation\")\n",
+ "p3 = plot_errors(steps, err_x1, h_star, \"cos(x1)\")\n",
+ "\n",
+ "plot(p1, p2, p3, layout=(1,3), size=(850, 350))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.9.0",
+ "language": "julia",
+ "name": "julia-1.9"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.9.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/tp/simple-shooting.ipynb b/tp/simple-shooting.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..95eb36385c58c711dcad707cc724d88a5088bdbc
--- /dev/null
+++ b/tp/simple-shooting.ipynb
@@ -0,0 +1,368 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "efcf7482",
+ "metadata": {},
+ "source": [
+ "[<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)\n",
+ "\n",
+ "# Tir simple indirect\n",
+ "\n",
+ "- Date : 2023\n",
+ "- Durée approximative : 3/4 de séance\n",
+ "\n",
+ "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.\n",
+ "\n",
+ "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",
+ "execution_count": null,
+ "id": "1e3057d1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Packages\n",
+ "using DifferentialEquations\n",
+ "using Plots\n",
+ "using NLsolve\n",
+ "using ForwardDiff\n",
+ "using LinearAlgebra\n",
+ "include(\"utils.jl\"); # plot_traj!, plot_flow!, Flow"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "1b19260f",
+ "metadata": {},
+ "source": [
+ "## Exercice 1\n",
+ "\n",
+ "On considère le problème de contrôle optimal suivant :\n",
+ "\n",
+ "$$ \n",
+ " \\left\\{ \n",
+ " \\begin{array}{l}\n",
+ " \\displaystyle \\min\\, \\frac{1}{2} \\int_0^{1} u^2(t) \\, \\mathrm{d}t \\\\[1.0em]\n",
+ " \\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]\n",
+ " x(0) = -1, \\quad x(1) = 0.\n",
+ " \\end{array}\n",
+ " \\right. \n",
+ "$$\n",
+ "\n",
+ "Pour $\\alpha \\in \\{0, 1.5\\}$, faire :\n",
+ "\n",
+ "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. \n",
+ "\n",
+ "**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`.\n",
+ "\n",
+ "2. Coder la fonction de tir \n",
+ "\n",
+ "$$\n",
+ "S(p_0) = \\pi(z(t_f, x_0, p_0)) - x_f,\n",
+ "$$\n",
+ "\n",
+ "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)`.\n",
+ "\n",
+ "3. Afficher la fonction de tir pour $p_0 \\in [0, 1]$.\n",
+ "\n",
+ "4. Résoudre l'équation de tir $S(p_0) = 0$ et tracer l'extrémale dans le plan de phase, cf. question 1.\n",
+ "\n",
+ "**Remarque.** Donner une bonne initialisation au solveur de Newton pour la résolution de l'équation de tir."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "47d35d5b-b659-4906-bd1a-61e369de876a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Paramètres du problème\n",
+ "α = 0\n",
+ "t0 = 0\n",
+ "tf = 1\n",
+ "x0 = -1\n",
+ "x_target = 0; # correspond à xf dans la description mathématique du problème"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "497d73ec-5b23-459e-8c32-ca13add12f5f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 1:\n",
+ "\n",
+ "# ------------------------------------\n",
+ "# conseils : \n",
+ "# - revoir la fin du TP précédent, ie le TP ode\n",
+ "# - réutiliser plot_traj! et plot_flow! comme au TP ode, cf utils.jl\n",
+ "#\n",
+ "\n",
+ "# ------------------------------------\n",
+ "# système hamiltonien sous la forme F(z), z = (x, p)\n",
+ "# \n",
+ "# à compléter ici\n",
+ "\n",
+ "# fin à compléter ici\n",
+ "# ------------------------------------\n",
+ "\n",
+ "# plot initial et paramètres du plot\n",
+ "plt = plot()\n",
+ "ymin = -0.2\n",
+ "ymax = 2.5\n",
+ "\n",
+ "# ------------------------------------\n",
+ "# plot de trajectoires\n",
+ "# \n",
+ "# à compléter ici\n",
+ "\n",
+ "# fin à compléter ici\n",
+ "# ------------------------------------\n",
+ "\n",
+ "\n",
+ "# ------------------------------------\n",
+ "# plot de fronts\n",
+ "# \n",
+ "# à compléter ici\n",
+ "\n",
+ "# fin à compléter ici\n",
+ "# ------------------------------------\n",
+ "\n",
+ "# additional plots\n",
+ "plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax), legend = false)\n",
+ "plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)\n",
+ "plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(900, 600), xlabel = \"x\", ylabel = \"p\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "df2c37e6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 2 : coder la fonction de tir\n",
+ "\n",
+ "# Système hamiltonien\n",
+ "function Hv(x, p)\n",
+ " \"\"\" \n",
+ " Calcul de la dynamique hamiltonienne\n",
+ " -------\n",
+ " \n",
+ " parameters (input)\n",
+ " ------------------\n",
+ " x : état\n",
+ " p : état adjoint\n",
+ " \n",
+ " returns\n",
+ " -------\n",
+ " dx = ẋ, dp = ṗ\n",
+ " \"\"\"\n",
+ " dx = 0\n",
+ " dp = 0\n",
+ " return dx, dp\n",
+ "end\n",
+ "\n",
+ "# flot du système hamiltonien: voir utils.jl\n",
+ "f = Flow(Hv) \n",
+ "\n",
+ "# fonction de tir\n",
+ "function shoot(p0)\n",
+ " return 0\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ba8b9913-f4f7-4a2f-9003-8179bf157aea",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 3 : afficher la fonction de tir\n",
+ "\n",
+ "# Plot initial\n",
+ "plt_2 = plot(size=(600, 400),xlabel = \"p₀\", ylabel = \"S(p₀)\", xticks = 0:0.1:1, xlims = (0,1), ylims = (-0.5,1) )\n",
+ "plot!(plt_2, [0,1], [0,0], legend = false, color = :black)\n",
+ "\n",
+ "p0 = range(0, 1, length = 50)\n",
+ "\n",
+ "#\n",
+ "# compléter à partir d'ici finaliser l'affichage\n",
+ "\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ea2ab41c-44d1-4e24-a99a-32f7e2a0b75a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 4 : résoudre S(p0) = 0\n",
+ "\n",
+ "# transformation de la fonction de tir en une fonction vectorielle: contrainte du solveur\n",
+ "S(p0) = [shoot(p0[1])]\n",
+ "\n",
+ "# jacobienne de la fonction de tir\n",
+ "JS(p0) = ForwardDiff.jacobian(S, p0)\n",
+ "\n",
+ "# initialisation de la fonction de tir \n",
+ "p0_guess = [NaN] # A MODIFIER\n",
+ "\n",
+ "# résolution de shoot(p0) = 0. L'itéré initial est appelé p0_guess\n",
+ "sol = nlsolve(S, JS, p0_guess; xtol=1e-8, method=:trust_region, show_trace=true);\n",
+ "println(\"\\n\", sol, \"\\n\")\n",
+ "\n",
+ "# récupération de la solution \n",
+ "p0_sol = sol.zero[1]\n",
+ "\n",
+ "#\n",
+ "println(\"p0 (solution) = \", p0_sol)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fe60edc8-df13-454e-a12d-700ec7f2419c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Affichage de la BC-extrémale\n",
+ "#\n",
+ "# CONSEIL : utiliser la fonction f, cf. f = Flow(Hv). Voir aussi utils.jl.\n",
+ "#\n",
+ "plt # A REMPLACER: attention, il faut remettre les xlims et ylims pour un affichage propre"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "15fe1891",
+ "metadata": {},
+ "source": [
+ "## Exercice 2\n",
+ "\n",
+ "On considère le problème de contrôle optimal suivant :\n",
+ "\n",
+ "$$ \n",
+ " \\left\\{ \n",
+ " \\begin{array}{l}\n",
+ " \\displaystyle \\min\\, \\frac{1}{2} \\int_0^{t_f} ||u(t)||^2 \\, \\mathrm{d}t \\\\[1.0em]\n",
+ " \\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]\n",
+ " x(0) = a, \\quad x(t_f) \\in b + \\mathrm{R} v,\n",
+ " \\end{array}\n",
+ " \\right. \n",
+ "$$\n",
+ "\n",
+ "avec $a = (0,0)$, $b = (1,0)$, $v = (0, 1)$ et $t_f = 1$.\n",
+ "\n",
+ "Pour différente valeur de $\\alpha \\in [0, 2]$, faire :\n",
+ "\n",
+ "1. Résoudre le problème de contrôle optimal par du tir simple.\n",
+ "\n",
+ "2. Tracer la trajectoire solution dans le plan $(x_1, x_2)$.\n",
+ "Tracer pour la solution, attaché au point $x(t_f)$, les vecteurs $\\dot{x}(t_f)$ et $p(t_f)$. Commentaires ?\n",
+ "\n",
+ "**Remarque.** Utiliser la fonction suivante pour tracer un vecteur $v = (v_1, v_2)$ aux coordonnées $(a, b)$ :\n",
+ "\n",
+ "``` julia \n",
+ " quiver!(plt, [a], [b], quiver=([v1], [v2]), c=\"green\", linewidth=2)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2a0f0b8b-5092-41cb-bb53-815c99ea95ee",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Paramètres du problème\n",
+ "α = 2\n",
+ "t0 = 0\n",
+ "tf = 1\n",
+ "x0 = [0, 0];"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "90248383",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 1 : résolution de S(p0) = 0\n",
+ "\n",
+ "# système hamiltonien \n",
+ "function Hv(x, p)\n",
+ " return [0, 0], [0, 0]\n",
+ "end\n",
+ "\n",
+ "# flot du système hamiltonien\n",
+ "f = Flow(Hv)\n",
+ "\n",
+ "# fonction de tir\n",
+ "function shoot(p0)\n",
+ " return [0, 0]\n",
+ "end\n",
+ "\n",
+ "# jacobienne de la fonction de tir\n",
+ "jshoot(p0) = ForwardDiff.jacobian(shoot, p0)\n",
+ "\n",
+ "# initialisation de la fonction de tir\n",
+ "p0_guess = [0.1, 0.1]\n",
+ "\n",
+ "# résolution de shoot(p0) = 0. L'itéré initial est appelé p0_guess\n",
+ "sol = nlsolve(shoot, jshoot, p0_guess; xtol=1e-8, method=:trust_region, show_trace=true);\n",
+ "println(\"\\n\", sol, \"\\n\")\n",
+ "\n",
+ "# récupération de la solution\n",
+ "p0_sol = sol.zero\n",
+ "\n",
+ "#\n",
+ "println(\"p0 (solution) = \", p0_sol)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8bd070d1-8474-4fbb-b0d5-054fdc5bffcd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "## Question 2 : affichage de la solution\n",
+ "\n",
+ "# calcul de la solution\n",
+ "ode_sol = f((t0, tf), x0, p0_sol);\n",
+ "\n",
+ "# affichage de la contrainte terminale\n",
+ "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\n",
+ "annotate!(plt, 1.05, -1, text(\"c=0\", :left, :blue, 12))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.8.2",
+ "language": "julia",
+ "name": "julia-1.8"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.8.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}