From aee635128cd032ea9f77f930310fade1f6b553bc Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Fri, 16 May 2025 14:42:20 +0200
Subject: [PATCH] up tp

---
 README.md                                     |   6 +-
 tp/derivative.ipynb                           | 495 ++++++++++++++++++
 tp/mpc-navigation.ipynb                       | 493 +++++++++++++++++
 tp/{runge-kutta.ipynb => rk-explicites.ipynb} |   0
 4 files changed, 992 insertions(+), 2 deletions(-)
 create mode 100644 tp/derivative.ipynb
 create mode 100644 tp/mpc-navigation.ipynb
 rename tp/{runge-kutta.ipynb => rk-explicites.ipynb} (100%)

diff --git a/README.md b/README.md
index dbf2dda..363606b 100644
--- a/README.md
+++ b/README.md
@@ -18,10 +18,12 @@ 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
+* [Sujet 4 : tp/mpc-navigation.ipynb](tp/mpc-navigation.ipynb) - Méthode MPC pour la navigation marine
 
-Equation différentielle ordinaire
+Calcul différentiel et équation différentielle ordinaire
 
-* [Sujet 4 : tp/runge-kutta.ipynb](tp/runge-kutta.ipynb) - Méthodes de Runge-Kutta
+* [Sujet 5 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul de dérivées
+* [Sujet 6 : tp/runge-kutta.ipynb](tp/runge-kutta.ipynb) - Méthodes de Runge-Kutta
 
 **Notes de cours supplémentaires - ressources**
 
diff --git a/tp/derivative.ipynb b/tp/derivative.ipynb
new file mode 100644
index 0000000..c03c237
--- /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 : 2025\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/mpc-navigation.ipynb b/tp/mpc-navigation.ipynb
new file mode 100644
index 0000000..9942438
--- /dev/null
+++ b/tp/mpc-navigation.ipynb
@@ -0,0 +1,493 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "660275fc",
+   "metadata": {},
+   "source": [
+    "<div style=\"display:flex\"> \n",
+    "    <img \n",
+    "      src=\"https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png\" \n",
+    "      height=\"100\"\n",
+    "    > \n",
+    "    <img \n",
+    "      src=\"https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/ship.jpg\"\n",
+    "      style=\"display: block;\n",
+    "      margin-left: auto;\n",
+    "      margin-right: auto;\n",
+    "      width: 50%;\"\n",
+    "    >\n",
+    "</div> \n",
+    "\n",
+    "# Problème de navigation, une approche [MPC](https://en.wikipedia.org/wiki/Model_predictive_control)\n",
+    "\n",
+    "- Date : 2025\n",
+    "- Durée approximative : inconnue\n",
+    "\n",
+    "On considère un navire dans un courant constant $w=(w_x,w_y)$, $\\|w\\| \\lt 1$. \n",
+    "[L'angle de cap](https://fr.wikipedia.org/wiki/Cap_(navigation)) est contrôlé, amenant aux équations différentielles suivantes : \n",
+    "\n",
+    "$$ \\begin{array}{rcl}\n",
+    "     \\dot{x}(t) &=& w_x+\\cos\\theta(t),\\quad t \\in [0,t_f]\\\\\n",
+    "     \\dot{y}(t) &=& w_y+\\sin\\theta(t),\\\\\n",
+    "     \\dot{\\theta}(t) &=& u(t). \n",
+    "   \\end{array} $$\n",
+    "\n",
+    "La vitesse angulaire est limitée et normalisée : $\\|u(t)\\| \\leq 1$. Il y a des conditions aux limites au temps initial $t = 0$ et au temps final $t = t_f$, sur la position $(x,y)$ et sur l'angle $\\theta$. L'objectif est de minimiser le temps final. Ce sujet est inspiré de ce [TP](https://github.com/pns-mam/commande/tree/main/tp3) dont le problème vient d'une collaboration entre l'Université Côte d'Azur et l'entreprise française [CGG](https://www.cgg.com) qui s'intéresse aux manoeuvres optimales de très gros navires pour la prospection marine.\n",
+    "\n",
+    "✏️ **Exercice 1.** On supposera $p^0 = -1$, on se place dans le cas normal.\n",
+    "\n",
+    "1. Ecrire le problème de contrôle optimal sous la forme de Mayer.\n",
+    "2. Donner le pseudo-hamiltonien $H(q, p, u)$, où $q = (x, y, \\theta)$ et $p = (p_x, p_y, p_\\theta)$.\n",
+    "3. Calculer l'équation adjointe, c'est-à-dire vérifiée par le vecteur adjoint $p$, donnée par le principe du maximum de Pontryagin.\n",
+    "4. Calculer le contrôle maximisant en fonction de $p_\\theta$ (on pourra l'écrire comme une [fonction multivaluée](https://fr.wikipedia.org/wiki/Fonction_multivaluée)).\n",
+    "5. Calculer le contrôle singulier, c'est-à-dire celui permettant de vérifier $p_\\theta(t) = 0$ sur un intervalle de temps non réduit à un singleton."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a68c170f",
+   "metadata": {},
+   "source": [
+    "## Données du problème"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8dd6f718",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OptimalControl, NLPModelsIpopt, Plots, OrdinaryDiffEq, LinearAlgebra, Plots.PlotMeasures\n",
+    "\n",
+    "t0 = 0.\n",
+    "x0 = 0. \n",
+    "y0 = 0.\n",
+    "θ0 = π/7\n",
+    "xf = 4.\n",
+    "yf = 7.\n",
+    "θf = -π/2\n",
+    "\n",
+    "function current(x, y) # current as a function of position\n",
+    "    ε = 1e-1\n",
+    "    w = [ 0.6, 0.4 ]\n",
+    "    δw = ε * [ y, -x ] / sqrt(1+x^2+y^2)\n",
+    "    w = w + δw\n",
+    "    if (w[1]^2 + w[2]^2 >= 1)\n",
+    "        error(\"|w| >= 1\")\n",
+    "    end\n",
+    "    return w\n",
+    "end\n",
+    "\n",
+    "#\n",
+    "function plot_state!(plt, x, y, θ; color=1)\n",
+    "    plot!(plt, [x], [y], marker=:circle, legend=false, color=color, markerstrokecolor=color, markersize=5, z_order=:front)\n",
+    "    quiver!(plt, [x], [y], quiver=([cos(θ)], [sin(θ)]), color=color, linewidth=2, z_order=:front)\n",
+    "    return plt\n",
+    "end\n",
+    "\n",
+    "function plot_current!(plt; current=current, N=10, scaling=1)\n",
+    "    for x ∈ range(xlims(plt)..., N)\n",
+    "        for y ∈ range(ylims(plt)..., N)\n",
+    "            w = scaling*current(x, y)\n",
+    "            quiver!(plt, [x], [y], quiver=([w[1]], [w[2]]), color=:black, linewidth=0.5, z_order=:back)\n",
+    "        end\n",
+    "    end\n",
+    "    return plt\n",
+    "end\n",
+    "\n",
+    "# on affiche dans le plan de phase augmenté les conditions aux limites et le courant\n",
+    "plt = plot(\n",
+    "    xlims=(-2, 6), \n",
+    "    ylims=(-1, 8), \n",
+    "    size=(600, 600), \n",
+    "    aspect_ratio=1, \n",
+    "    xlabel=\"x\", \n",
+    "    ylabel=\"y\", \n",
+    "    title=\"Conditions aux limites\",\n",
+    "    leftmargin=5mm, \n",
+    "    bottommargin=5mm,\n",
+    ")\n",
+    "\n",
+    "plot_state!(plt, x0, y0, θ0; color=2)\n",
+    "plot_state!(plt, xf, yf, θf; color=2)\n",
+    "plot_current!(plt)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b4c8930",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "function plot_trajectory!(plt, t, x, y, θ; N=5) # N : nombre de points où l'on va afficher θ\n",
+    "\n",
+    "    # trajectoire\n",
+    "    plot!(plt, x.(t), y.(t), legend=false, color=1, linewidth=2, z_order=:front)\n",
+    "\n",
+    "    if N > 0\n",
+    "\n",
+    "        # longueur du trajet\n",
+    "        s = 0\n",
+    "        for i ∈ 2:length(t)\n",
+    "            s += norm([x(t[i]), y(t[i])] - [x(t[i-1]), y(t[i-1])])\n",
+    "        end\n",
+    "\n",
+    "        # intervalle de longueur\n",
+    "        Δs = s/(N+1)\n",
+    "        tis = []\n",
+    "        s = 0\n",
+    "        for i ∈ 2:length(t)\n",
+    "            s += norm([x(t[i]), y(t[i])] - [x(t[i-1]), y(t[i-1])])\n",
+    "            if s > Δs && length(tis) < N\n",
+    "                push!(tis, t[i])\n",
+    "                s = 0\n",
+    "            end\n",
+    "        end\n",
+    "\n",
+    "        # affichage des points intermédiaires\n",
+    "        for ti ∈ tis\n",
+    "            plot_state!(plt, x(ti), y(ti), θ(ti); color=1)\n",
+    "        end\n",
+    "\n",
+    "    end\n",
+    "\n",
+    "    return plt\n",
+    "    \n",
+    "end;"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0403331e",
+   "metadata": {},
+   "source": [
+    "## Solveur (OptimalControl)\n",
+    "\n",
+    "✏️ **Exercice 2.** Coder le problème de contrôle optimal ci-dessous.\n",
+    "\n",
+    "On pourra s'insiprer du problème de la [rame de métro à temps minimal](https://control-toolbox.org/OptimalControl.jl/stable/tutorial-double-integrator-time.html) décrit dans la documentation du package OptimalControl pour définir notre problème de manoeuvre de navire ci-après.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1933e2c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "function solve(t0, x0, y0, θ0, xf, yf, θf, w; \n",
+    "    grid_size=300, tol=1e-8, max_iter=500, print_level=4, display=true)\n",
+    "\n",
+    "    # ---------------------------------------\n",
+    "    # Définition du problème : TO UPDATE\n",
+    "    ocp = @def begin\n",
+    "\n",
+    "    end\n",
+    "    # ---------------------------------------\n",
+    "\n",
+    "    # Initialisation\n",
+    "    tf_init = 1.5*norm([xf, yf]-[x0, y0])\n",
+    "    x_init(t) = [ x0, y0, θ0 ] * (tf_init-t)/(tf_init-t0) + [xf, yf, θf] * (t-t0)/(tf_init-t0)\n",
+    "    u_init = (θf - θ0) / (tf_init-t0)\n",
+    "    init = (state=x_init, control=u_init, variable=tf_init)\n",
+    "\n",
+    "    # Résolution\n",
+    "    sol = OptimalControl.solve(ocp; \n",
+    "        init=init,\n",
+    "        grid_size=grid_size, \n",
+    "        tol=tol, \n",
+    "        max_iter=max_iter, \n",
+    "        print_level=print_level,\n",
+    "        display=display,\n",
+    "        disc_method=:euler,\n",
+    "    )\n",
+    "\n",
+    "    # Récupération des données utiles\n",
+    "    t = time_grid(sol)\n",
+    "    q = state(sol)\n",
+    "    x = t -> q(t)[1]\n",
+    "    y = t -> q(t)[2]\n",
+    "    θ = t -> q(t)[3]\n",
+    "    u = control(sol)\n",
+    "    tf = variable(sol)\n",
+    "    \n",
+    "    return t, x, y, θ, u, tf, iterations(sol), sol.solver_infos.constraints_violation\n",
+    "    \n",
+    "end;"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cc8547d2",
+   "metadata": {},
+   "source": [
+    "## Première résolution\n",
+    "\n",
+    "On considère ici un vent constant et on résout une première fois le problème."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ffea082a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# -----------------------------------------------------------------------------------------------------------\n",
+    "t, x, y, θ, u, tf, iter, cons = solve(t0, x0, y0, θ0, xf, yf, θf, current(x0, y0); display=false);\n",
+    "\n",
+    "println(\"Iterations : \", iter)\n",
+    "println(\"Constraints violation : \", cons)\n",
+    "println(\"tf : \", tf)\n",
+    "\n",
+    "# -----------------------------------------------------------------------------------------------------------\n",
+    "# affichage\n",
+    "\n",
+    "# trajectoire\n",
+    "plt_q = plot(xlims=(-2, 6), ylims=(-1, 8), aspect_ratio=1, xlabel=\"x\", ylabel=\"y\")\n",
+    "plot_state!(plt_q, x0, y0, θ0; color=2)\n",
+    "plot_state!(plt_q, xf, yf, θf; color=2)\n",
+    "plot_current!(plt_q; current=(x, y) -> current(x0, y0))\n",
+    "plot_trajectory!(plt_q, t, x, y, θ)\n",
+    "\n",
+    "# contrôle\n",
+    "plt_u = plot(t, u; color=1, legend=false, linewidth=2, xlabel=\"t\", ylabel=\"u\")\n",
+    "\n",
+    "#\n",
+    "plot(plt_q, plt_u; \n",
+    "    layout=(1, 2), \n",
+    "    size=(1200, 600),\n",
+    "    leftmargin=5mm, \n",
+    "    bottommargin=5mm,\n",
+    "    plot_title=\"Simulation courant constant\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "827b18d5",
+   "metadata": {},
+   "source": [
+    "## Simulation du système réel\n",
+    "\n",
+    "Dans la simulation précédente, nous faisons l'hypothèse que le courant est constant. Cependant, d'un point de vue pratique le courant dépend de la position $(x, y)$. Etant donné un modèle de courant, donné par la fonction `current`, nous pouvons simuler la trajectoire réelle du navire, pourvu que l'on ait la condition initiale et le contrôle au cours du temps."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "25422e2e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "function realistic_trajectory(tf, t0, x0, y0, θ0, u, current; abstol=1e-12, reltol=1e-12, saveat=[])\n",
+    "    \n",
+    "    function rhs!(dq, q, dummy, t)\n",
+    "        x, y, θ = q\n",
+    "        w = current(x, y)\n",
+    "        dq[1] = w[1]+cos(θ)\n",
+    "        dq[2] = w[2]+sin(θ)\n",
+    "        dq[3] = u(t)\n",
+    "    end\n",
+    "    \n",
+    "    q0 = [ x0, y0, θ0 ]\n",
+    "    tspan = (t0, tf)\n",
+    "    ode = ODEProblem(rhs!, q0, tspan)\n",
+    "    sol = OrdinaryDiffEq.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)\n",
+    "\n",
+    "    t = sol.t\n",
+    "    x = t -> sol(t)[1]\n",
+    "    y = t -> sol(t)[2]\n",
+    "    θ = t -> sol(t)[3]\n",
+    "\n",
+    "    return t, x, y, θ\n",
+    "    \n",
+    "end;"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da61110e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# trajectoire réaliste\n",
+    "t, x, y, θ = realistic_trajectory(tf, t0, x0, y0, θ0, u, current)\n",
+    "\n",
+    "# -----------------------------------------------------------------------------------------------------------\n",
+    "# affichage\n",
+    "\n",
+    "# trajectoire\n",
+    "plt_q = plot(xlims=(-2, 6), ylims=(-1, 8), aspect_ratio=1, xlabel=\"x\", ylabel=\"y\")\n",
+    "plot_state!(plt_q, x0, y0, θ0; color=2)\n",
+    "plot_state!(plt_q, xf, yf, θf; color=2)\n",
+    "plot_current!(plt_q; current=current)\n",
+    "plot_trajectory!(plt_q, t, x, y, θ)\n",
+    "plot_state!(plt_q, x(tf), y(tf), θ(tf); color=3)\n",
+    "\n",
+    "# contrôle\n",
+    "plt_u = plot(t, u; color=1, legend=false, linewidth=2, xlabel=\"t\", ylabel=\"u\")\n",
+    "\n",
+    "#\n",
+    "plot(plt_q, plt_u; \n",
+    "    layout=(1, 2), \n",
+    "    size=(1200, 600),\n",
+    "    leftmargin=5mm, \n",
+    "    bottommargin=5mm,\n",
+    "    plot_title=\"Simulation avec modèle de courant\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "503bd360",
+   "metadata": {},
+   "source": [
+    "## Approche MPC\n",
+    "\n",
+    "En pratique, nous n'avons pas à l'avance les données réelles du courant sur l'ensemble du trajet, c'est pourquoi nous allons recalculer régulièrement le contrôle optimal. L'idée est de mettre à jour le contrôle optimal à intervalle de temps régulier en prenant en compte le courant à la position où le navire se trouve. On est donc amener à résoudre un certain nombre de problème à courant constant, avec celui-ci mis réguilièremet à jour. Ceci est une introduction aux méthodes dites MPC, pour \"Model Predictive Control\" en anglais."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "563c73fd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Nmax = 20   # nombre d'itérations max de la méthode MPC\n",
+    "ε = 1e-1    # rayon sur la condition finale pour stopper les calculs\n",
+    "Δt = 1.0    # pas de temps fixe de la méthode MPC\n",
+    "P = 300      # nombre de points de discrétisation pour le solveur\n",
+    "\n",
+    "t1 = t0\n",
+    "x1 = x0\n",
+    "y1 = y0\n",
+    "θ1 = θ0\n",
+    "\n",
+    "data = []\n",
+    "\n",
+    "N = 1\n",
+    "stop = false\n",
+    "\n",
+    "while !stop\n",
+    "    \n",
+    "    # on récupère le courant à la position actuelle\n",
+    "    w = current(x1, y1)\n",
+    "\n",
+    "    # on résout le problème\n",
+    "    t, x, y, θ, u, tf, iter, cons = solve(t1, x1, y1, θ1, xf, yf, θf, w; grid_size=P, display=false);\n",
+    "\n",
+    "    # calcul du temps suivant\n",
+    "    if (t1+Δt < tf)\n",
+    "        t2 = t1+Δt\n",
+    "    else\n",
+    "        t2 = tf\n",
+    "        println(\"t2=tf: \", t2)\n",
+    "        stop = true\n",
+    "    end\n",
+    "\n",
+    "    # on stocke les données: le temps initial courant, le temps suivant, le contrôle\n",
+    "    push!(data, (t2, t1, x(t1), y(t1), θ(t1), u, tf))\n",
+    "\n",
+    "    # on met à jour les paramètres de la méthode MPC: on simule la réalité\n",
+    "    t, x, y, θ = realistic_trajectory(t2, t1, x1, y1, θ1, u, current)\n",
+    "    t1 = t2\n",
+    "    x1 = x(t1)\n",
+    "    y1 = y(t1)\n",
+    "    θ1 = θ(t1)\n",
+    "\n",
+    "    # on calcule la distance à la position cible\n",
+    "    distance = norm([ x1, y1, θ1 ] - [ xf, yf, θf ])\n",
+    "    println(\"N: \", N, \"\\t distance: \", distance, \"\\t iterations: \", iter, \"\\t constraints: \", cons, \"\\t tf: \", tf)\n",
+    "    if !((distance > ε) && (N < Nmax))\n",
+    "        stop = true\n",
+    "    end\n",
+    "\n",
+    "    #\n",
+    "    N += 1\n",
+    "\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "edd803b2",
+   "metadata": {},
+   "source": [
+    "## Affichage"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d0292a49",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# trajectoire\n",
+    "plt_q = plot(xlims=(-2, 6), ylims=(-1, 8), aspect_ratio=1, xlabel=\"x\", ylabel=\"y\")\n",
+    "\n",
+    "# condition finale\n",
+    "plot_state!(plt_q, xf, yf, θf; color=2)\n",
+    "\n",
+    "# courant\n",
+    "plot_current!(plt_q; current=current)\n",
+    "\n",
+    "# contrôle\n",
+    "plt_u = plot(xlabel=\"t\", ylabel=\"u\")\n",
+    "\n",
+    "N = 1\n",
+    "\n",
+    "for d ∈ data\n",
+    "\n",
+    "    #\n",
+    "    t2, t1, x1, y1, θ1, u, tf = d\n",
+    "\n",
+    "    # calcule de la trajectoire réelle\n",
+    "    t, x, y, θ = realistic_trajectory(t2, t1, x1, y1, θ1, u, current)\n",
+    "\n",
+    "    # trajectoire\n",
+    "    plot_state!(plt_q, x1, y1, θ1; color=2)\n",
+    "    plot_trajectory!(plt_q, t, x, y, θ; N=0)\n",
+    "\n",
+    "    # contrôle\n",
+    "    plot!(plt_u, t, u; color=1, legend=false, linewidth=2)\n",
+    "\n",
+    "    N += 1\n",
+    "\n",
+    "end\n",
+    "\n",
+    "plot_state!(plt_q, x(tf), y(tf), θ(tf); color=3)\n",
+    "\n",
+    "#\n",
+    "plot(plt_q, plt_u; \n",
+    "    layout=(1, 2), \n",
+    "    size=(1200, 600),\n",
+    "    leftmargin=5mm, \n",
+    "    bottommargin=5mm,\n",
+    "    plot_title=\"Simulation avec modèle de courant\"\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.11.3",
+   "language": "julia",
+   "name": "julia-1.11"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.11.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tp/runge-kutta.ipynb b/tp/rk-explicites.ipynb
similarity index 100%
rename from tp/runge-kutta.ipynb
rename to tp/rk-explicites.ipynb
-- 
GitLab