diff --git a/README.md b/README.md
index 0d629a560729c63108f7aada53119dbd56bf97d1..7ef0c340f81e1baffbf5a2c307e39fff235f0681 100644
--- a/README.md
+++ b/README.md
@@ -14,6 +14,16 @@ 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)
diff --git a/tp/derivative.ipynb b/tp/derivative.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..c03c237a6c51ebd38d388e582fb8c73022548d2f
--- /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/direct-indirect.ipynb b/tp/direct-indirect.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..003ebfa603ea437320f7a8f01e067ae7bf45f29c
--- /dev/null
+++ b/tp/direct-indirect.ipynb
@@ -0,0 +1,459 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "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",
+ "# Méthode directe et indirecte\n",
+ "\n",
+ "- Date : 2025\n",
+ "- Durée approximative : 1.5 séance\n",
+ "\n",
+ "Le but est de résoudre par une méthode de tir, 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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "using JuMP, Ipopt # pour la méthode directe\n",
+ "using DifferentialEquations, NLsolve, ForwardDiff # pour la méthode indirecte\n",
+ "using Plots # pour les graphiques\n",
+ "include(\"utils.jl\"); # fonctions utilitaires"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<div style=\"width:95%;\n",
+ " margin:10px;\n",
+ " padding:8px;\n",
+ " color:white;\n",
+ " background-color: rgb(46, 109, 4);\n",
+ " border-radius:10px;\n",
+ " font-weight:bold;\n",
+ " font-size:1.5em;\n",
+ " text-align:center;\">\n",
+ "Introduction\n",
+ "</div>\n",
+ "\n",
+ "On considère le problème de contrôle optimal suivant :\n",
+ "\n",
+ "$$ \n",
+ " \\left\\{ \n",
+ " \\begin{array}{l}\n",
+ " \\displaystyle \\min_{x, u} \\displaystyle \\int_0^{t_f} x^2(t) \\, \\mathrm{d}t \\\\[1.0em]\n",
+ " \\dot{x}(t) = \\displaystyle u(t), \\quad |u(t)| \\le 1, \\quad t \\in [0, t_f] \\text{ p.p.}, \\\\[1.0em]\n",
+ " x(0) = 1, \\quad x(t_f) = 1/2\n",
+ " \\end{array}\n",
+ " \\right. \n",
+ "$$\n",
+ "où $t_f$ = 2.\n",
+ "\n",
+ "✏️ **Exercice 1.**\n",
+ "\n",
+ "1. Appliquez le PMP. Que pouvez-vous dire du contrôle maximisant ? \n",
+ "2. Peut-on appliquer simplement la méthode de tir simple vu au TP précédent ?"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<div style=\"width:95%;\n",
+ " margin:10px;\n",
+ " padding:8px;\n",
+ " color:white;\n",
+ " background-color: rgb(46, 109, 4);\n",
+ " border-radius:10px;\n",
+ " font-weight:bold;\n",
+ " font-size:1.5em;\n",
+ " text-align:center;\">\n",
+ "Méthode directe\n",
+ "</div>\n",
+ "\n",
+ "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.\n",
+ "\n",
+ "✏️ **Exercice 2.** \n",
+ "\n",
+ "- Mettez le problème sous forme de Mayer (cf. cours). Vous nommerez la nouvelle variable d'état associée au coût $c(\\cdot)$."
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Description de la méthode directe.**\n",
+ "\n",
+ "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. \n",
+ "\n",
+ "Pour cela :\n",
+ "\n",
+ "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.\n",
+ "2. On discrétise l'état, le contrôle et le coût sur cette grille et on note \n",
+ "$$\n",
+ "X = \\{ (x_i, u_i, c_i) ~|~ i \\in \\{1, \\dots, N+1\\}\\}\n",
+ "$$ \n",
+ "l'ensemble des variables d'optimisation du problème discrétisé.\n",
+ "\n",
+ "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\n",
+ "$$\n",
+ "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\\}.\n",
+ "$$\n",
+ "\n",
+ "✏️ ️**Exercice 3.** Définissez pour le problème discrétisé : \n",
+ "\n",
+ "1. l'objectif à minimiser en fonction d'une ou plusieurs composantes de $X$ bien choisies.\n",
+ "2. les contraintes d'inégalités associées au contrôle.\n",
+ "3. les contraintes initiales et finales associées à la variable d'état.\n",
+ "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)).\n"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "✏️ ️**Exercice 4.** \n",
+ "\n",
+ "- Modifiez le code suivant afin de résoudre le problème par la méthode directe que vous venez de décrire.\n",
+ "\n",
+ "**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",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Création du modèle JuMP, Utilisation de Ipopt comme solver\n",
+ "sys = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 5))\n",
+ "set_optimizer_attribute(sys,\"tol\",1e-8)\n",
+ "set_optimizer_attribute(sys,\"constr_viol_tol\",1e-6)\n",
+ "set_optimizer_attribute(sys,\"max_iter\",1000)\n",
+ "\n",
+ "#######\n",
+ "####### DEBUT A MODIFIER\n",
+ "#######\n",
+ "####### Les ... sont à remplacer !\n",
+ "#######\n",
+ "####### Attention, on écrit le problème sous la forme d'un problème de Mayer\n",
+ "#######\n",
+ "\n",
+ "# Paramètres\n",
+ "t0 = 0 # temps initial\n",
+ "tf = 0 # temps final\n",
+ "c0 = 0 # coût initial\n",
+ "x0 = 0 # état initial\n",
+ "xf = 0 # état final \n",
+ "\n",
+ "N = 50 # taille de la grille\n",
+ "Δt = (tf-t0)/N # pas de temps\n",
+ "\n",
+ "@variables(sys, begin\n",
+ " ... # coût\n",
+ " ... # état\n",
+ " ... # control\n",
+ "end)\n",
+ "\n",
+ "# Objectif\n",
+ "@objective(sys, Min, ...)\n",
+ "\n",
+ "# Conditions initiales et finales\n",
+ "@constraints(sys, begin\n",
+ " con_c0, ... # contraint sur le coût initial\n",
+ " con_x0, ... # contraint sur l'état initial\n",
+ " con_xf, ... # contraint sur l'état final\n",
+ "end)\n",
+ "\n",
+ "# Contraintes de dynamique avec le schéma de Crank-Nicolson\n",
+ "@NLconstraints(sys, begin\n",
+ " con_dc[j=1:N], ... # contraint différentielle sur le coût\n",
+ " con_dx[j=1:N], ... # contraint différentielle sur l'état\n",
+ "end);\n",
+ "\n",
+ "#######\n",
+ "####### FIN A MODIFIER\n",
+ "#######"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Solve for the control and state\n",
+ "println(\"Solving...\")\n",
+ "optimize!(sys)\n",
+ "println()\n",
+ "\n",
+ "# Display results\n",
+ "if termination_status(sys) == MOI.OPTIMAL\n",
+ " println(\" Solution is optimal\")\n",
+ "elseif termination_status(sys) == MOI.LOCALLY_SOLVED\n",
+ " println(\" (Local) solution found\")\n",
+ "elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys)\n",
+ " println(\" Solution is suboptimal due to a time limit, but a primal solution is available\")\n",
+ "else\n",
+ " error(\" The model was not solved correctly.\")\n",
+ "end\n",
+ "println(\" objective value = \", objective_value(sys))\n",
+ "println()\n",
+ "\n",
+ "# Retrieves values (including duals)\n",
+ "c = value.(c)[:]\n",
+ "x = value.(x)[:]\n",
+ "u = value.(u)[:]\n",
+ "t = (0:N) * value.(Δt)\n",
+ "\n",
+ "pc0 = dual(con_c0)\n",
+ "px0 = dual(con_x0)\n",
+ "pxf = dual(con_xf)\n",
+ "\n",
+ "if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end\n",
+ "if(px0*dual(con_dx[1])<0); px0 = -px0; end\n",
+ "if(pxf*dual(con_dx[N])<0); pxf = -pxf; end\n",
+ "\n",
+ "if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle\n",
+ " sign = -1.0\n",
+ "else\n",
+ " sign = 1.0\n",
+ "end\n",
+ "\n",
+ "pc = [ dual(con_dc[i]) for i in 1:N ]\n",
+ "px = [ dual(con_dx[i]) for i in 1:N ]\n",
+ "\n",
+ "pc = sign * [pc0; pc[1:N]] \n",
+ "px = sign * [px0; (px[1:N-1]+px[2:N])/2; pxf]; # We add the multiplier from the limit conditions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x_plot = plot(t, x, xlabel = \"t\", ylabel = \"x\", legend = false)\n",
+ "u_plot = plot(t, u, xlabel = \"t\", ylabel = \"u\", legend = false, size=(800,300), linetype=:steppre)\n",
+ "px_plot = plot(t, px, xlabel = \"t\", ylabel = \"p\", legend = false)\n",
+ "display(plot(x_plot, px_plot, layout = (1,2), size=(800,300)))\n",
+ "display(u_plot)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "✏️ ️**Exercice 5.** \n",
+ "\n",
+ "1. Commentez les résultats. \n",
+ "2. Modifiez la tolérance de l'optimiseur ainsi que le nombre de points de discrétisation. Commentaires.\n",
+ "3. Décrivez la structure du contrôle optimal en fonction du temps."
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<div style=\"width:95%;\n",
+ " margin:10px;\n",
+ " padding:8px;\n",
+ " color:white;\n",
+ " background-color: rgb(46, 109, 4);\n",
+ " border-radius:10px;\n",
+ " font-weight:bold;\n",
+ " font-size:1.5em;\n",
+ " text-align:center;\">\n",
+ "Méthode de tir indirect\n",
+ "</div>\n",
+ "\n",
+ "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.\n",
+ "\n",
+ "✏️ ️️**Exercice 6.** \n",
+ "\n",
+ "1. Donner la fonction de commutation.\n",
+ "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'est-à-dire le long d'un arc où la fonction de commutation est constante égale à 0.\n",
+ "3. Remplir le code ci-dessous: il faut fournir les 3 contrôles pour définir les flots, puis écrire la fonction de tir.\n",
+ "\n",
+ "**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.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# système pseudo-hamiltonien\n",
+ "function hv(x, p, u)\n",
+ " return [u, 2x]\n",
+ "end\n",
+ "\n",
+ "# systèmes hamiltoniens\n",
+ "hv_min(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur\n",
+ "hv_max(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur\n",
+ "hv_sin(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur\n",
+ "\n",
+ "# flots\n",
+ "fmin = Flow(hv_min)\n",
+ "fmax = Flow(hv_max)\n",
+ "fsin = Flow(hv_sin)\n",
+ "\n",
+ "# fonction de tir\n",
+ "function shoot(p0, t1, t2)\n",
+ " \n",
+ " # integration\n",
+ " x1, p1 = fmin(t0, x0, p0, t1) # x1, p1 sont des scalaires\n",
+ " x2, p2 = fsin(t1, x1, p1, t2)\n",
+ " x3, p3 = fmax(t2, x2, p2, tf)\n",
+ " \n",
+ " # conditions\n",
+ " s = zeros(eltype(p0), 3)\n",
+ " s[1] = NaN ## REMPLACER NaN par la bonne expression\n",
+ " s[2] = NaN ## REMPLACER NaN par la bonne expression\n",
+ " s[3] = NaN ## REMPLACER NaN par la bonne expression\n",
+ " \n",
+ " return s\n",
+ "\n",
+ "end;"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Itéré initial:\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[NaN, NaN, NaN]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# itéré initiale pour la méthode indirecte de tir\n",
+ "p0 = NaN ## REMPLACER NaN par une bonne valeur\n",
+ "t1 = NaN ## REMPLACER NaN par une bonne valeur\n",
+ "t2 = NaN ## REMPLACER NaN par une bonne valeur\n",
+ "\n",
+ "#\n",
+ "y = [ p0 ; t1 ; t2]\n",
+ "\n",
+ "println(\"Itéré initial:\\n\", y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# fonction de tir et sa jacobienne\n",
+ "foo(y) = shoot(y...)\n",
+ "jfoo(y) = ForwardDiff.jacobian(foo, y)\n",
+ "\n",
+ "# Résolution de shoot(p0, t1, t2) = 0.\n",
+ "nl_sol = nlsolve(foo, jfoo, y; xtol=1e-8, method=:trust_region, show_trace=true);\n",
+ "\n",
+ "# Retrieves solution\n",
+ "if converged(nl_sol)\n",
+ " p0 = nl_sol.zero[1]\n",
+ " t1 = nl_sol.zero[2]\n",
+ " t2 = nl_sol.zero[3];\n",
+ " println(\"\\nFinal solution:\\n\", nl_sol.zero)\n",
+ "else\n",
+ " error(\"Not converged\")\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# affichage de la solution\n",
+ "ode_sol = fmin((t0, t1), x0, p0)\n",
+ "tt0 = ode_sol.t\n",
+ "xx0 = [ ode_sol[1, j] for j in 1:size(tt0, 1) ]\n",
+ "pp0 = [ ode_sol[2, j] for j in 1:size(tt0, 1) ]\n",
+ "uu0 = -ones(size(tt0, 1))\n",
+ "\n",
+ "ode_sol = fsin((t1, t2), xx0[end], pp0[end])\n",
+ "tt1 = ode_sol.t\n",
+ "xx1 = [ ode_sol[1, j] for j in 1:size(tt1, 1) ]\n",
+ "pp1 = [ ode_sol[2, j] for j in 1:size(tt1, 1) ]\n",
+ "uu1 = zeros(size(tt1, 1))\n",
+ "\n",
+ "ode_sol = fmax((t2, tf), xx1[end], pp1[end])\n",
+ "tt2 = ode_sol.t\n",
+ "xx2 = [ ode_sol[1, j] for j in 1:size(tt2, 1) ]\n",
+ "pp2 = [ ode_sol[2, j] for j in 1:size(tt2, 1) ]\n",
+ "uu2 = ones(size(tt2, 1))\n",
+ "\n",
+ "t_shoot = [ tt0 ; tt1 ; tt2 ]\n",
+ "x_shoot = [ xx0 ; xx1 ; xx2 ]\n",
+ "p_shoot = [ pp0 ; pp1 ; pp2 ]\n",
+ "u_shoot = [ uu0 ; uu1 ; uu2 ] \n",
+ "\n",
+ "x_plot = plot(t_shoot, x_shoot, xlabel = \"t\", ylabel = \"x\", legend = false)\n",
+ "u_plot = plot(t_shoot, u_shoot, xlabel = \"t\", ylabel = \"u\", legend = false, size=(800,300), linetype=:steppre)\n",
+ "px_plot = plot(t_shoot, p_shoot, xlabel = \"t\", ylabel = \"px\", legend = false)\n",
+ "display(plot(x_plot, px_plot, layout = (1,2), size=(800, 300)))\n",
+ "display(u_plot)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "5. Comparer avec l'affichage de la solution de la méthode directe."
+ ]
+ }
+ ],
+ "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": 4
+}
diff --git a/tp/simple-shooting.ipynb b/tp/simple-shooting.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..c30f95da178520902b76531b80311afc22ee1b09
--- /dev/null
+++ b/tp/simple-shooting.ipynb
@@ -0,0 +1,367 @@
+{
+ "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 : 2025\n",
+ "- Durée approximative : 1.5 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",
+ "# 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 pour 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
+}