diff --git a/README.md b/README.md
index 7ef0c340f81e1baffbf5a2c307e39fff235f0681..0d629a560729c63108f7aada53119dbd56bf97d1 100644
--- a/README.md
+++ b/README.md
@@ -14,16 +14,6 @@ Systèmes dynamiques
* [Sujet 1 : tp/ode.ipynb](tp/ode.ipynb)
-Méthodes indirectes et directes en contrôle optimal
-
-* [Sujet 2 : tp/simple-shooting.ipynb](tp/simple-shooting.ipynb) - Tir simple indirect
-* [Sujet 3 : tp/direct-indirect.ipynb](tp/direct-indirect.ipynb) - Méthode directe et tir avec structure sur le contrôle
-
-Calcul différentiel et équation différentielle ordinaire
-
-* [Sujet 4 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul numérique de dérivées
-* [Sujet 5 : tp/runge-kutta.ipynb](tp/runge-kutta.ipynb) - Méthodes de Runge-Kutta
-
**Notes de cours supplémentaires - ressources**
* [Automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf)
diff --git a/tp/derivative.ipynb b/tp/derivative.ipynb
deleted file mode 100644
index f4869bb4a8c96caf4fc5c792afa7fa5da3d52cf2..0000000000000000000000000000000000000000
--- a/tp/derivative.ipynb
+++ /dev/null
@@ -1,495 +0,0 @@
-{
- "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/direct-indirect.ipynb b/tp/direct-indirect.ipynb
deleted file mode 100644
index 3e6a38e0f3c62501307d33d98168714d39387e55..0000000000000000000000000000000000000000
--- a/tp/direct-indirect.ipynb
+++ /dev/null
@@ -1,430 +0,0 @@
-{
- "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 tir simple structuré\n",
- "\n",
- "- Date : 2024\n",
- "- Durée approximative : 1.5 séance\n",
- "\n",
- "Le but est de résoudre par du tir simple indirect, un problème de contrôle optimal dont le contrôle est discontinu. On se propose de résoudre dans un premier temps, le problème par une méthode directe, afin de déterminer la structure optimale et une bonne approximation de la solution pour faire converger la méthode de tir."
- ]
- },
- {
- "cell_type": "code",
- "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": [
- "## Introduction\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": [
- "## Méthode directe\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 (c.f. 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",
- "#######\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": [
- "## Méthode indirecte\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-à-d 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/rk-implicites.ipynb b/tp/rk-implicites.ipynb
deleted file mode 100644
index 565a76874b1ef9f77fbc2e1631be59111717630f..0000000000000000000000000000000000000000
--- a/tp/rk-implicites.ipynb
+++ /dev/null
@@ -1,467 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png\" alt=\"N7\" height=\"80\"/>\n",
- "\n",
- "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"80\"/>\n",
- "\n",
- "# Méthodes de Runge-Kutta implicites - TP Projet\n",
- "\n",
- "- Date : 2023-2024\n",
- "- Durée approximative : inconnue"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "**Nom** : \"Entrez votre nom ici\"\n",
- "\n",
- "**Prénom** : \"Entrez votre prénom ici\""
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Rendu et consignes\n",
- "\n",
- "Une fois le travail terminé, vous enverrez directement le fichier `.ipynb` par email à l'adresse : `olivier.cots@toulouse-inp.fr`.\n",
- "\n",
- "- **Date limite du rendu : mercredi 15/11/2023 à 23h59.** Attention, à partir de 24h, 2 points est enlevé de la note finale toutes les 30 minutes.\n",
- "- **Attention :** Le fichier doit être renommé de la façon suivante : `rk_implicites_NOM_Prenom.ipynb`. 4 points enlevés si le nom du fichier n'est pas respecté.\n",
- "- **Documents autorisés :** vous pouvez utiliser votre polycopié et les TPs précédents."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "(-2, -2)"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "a = (1, 2)\n",
- "b = (3, 4)\n",
- "a.-b"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 17,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "f (generic function with 1 method)"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "f(x, y) = x+y "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 19,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "3"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "f(a...)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 16,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "2-element Vector{Int64}:\n",
- " 1\n",
- " 2"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "v = [a...]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 23,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "0.0"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "using LinearAlgebra\n",
- "k1 = [1, 2]\n",
- "k2 = [3, 4]\n",
- "ks = (k1, k2)\n",
- "G(k1, k2) = (k1, k2)\n",
- "norm(ks.-G(ks...))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Introduction\n",
- "\n",
- "Nous allons dans ce TP, implémenter quelques méthodes de Runge-Kutta **implicites** (voir **polycopié Section 8.2**) et étudier leur convergence. On considère un pas de temps $h$ uniforme. Une méthode à un pas implicite est convergente si pour toute solution $x(\\cdot, t_0, x_0)$ d'un problème de Cauchy, la suite approximante ${(x_i)}_i$ donnée par la méthode à un pas implicite vérifie \n",
- "$$\n",
- " \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\to 0 \n",
- " \\quad\\text{quand}\\quad h \\to 0.\n",
- "$$\n",
- "\n",
- "Si la convergence est d'ordre $p$ alors il existe une constante $C \\ge 0$ telle que l'**erreur globale** $E$ vérifie\n",
- "\n",
- "$$\n",
- " E := \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\le C\\, h^p.\n",
- "$$\n",
- "\n",
- "Faisons l'hypothèse que $E = M\\, h^p$ pour un certain $M \\ge 0$. En passant au logarithme, on obtient\n",
- "\n",
- "$$\n",
- " \\log(E) = \\log(M) + p\\, \\log(h).\n",
- "$$\n",
- "\n",
- "Nous en déduisons que si on trace $\\log(E)$ en fonction de $\\log(h)$, on doit obtenir une droite de pente $p$. C'est ce que nous allons vérifier dans ce TP."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# activate local project\n",
- "using Pkg\n",
- "Pkg.activate(\".\")\n",
- "\n",
- "# load packages\n",
- "using LinearAlgebra\n",
- "using Plots\n",
- "using Plots.PlotMeasures\n",
- "using Polynomials\n",
- "\n",
- "#\n",
- "px = PlotMeasures.px;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Fonctons auxiliaires \n",
- "\n",
- "# VOUS POUVEZ METTRE VOS FONCTIONS AUXILIAIRES ICI"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## L'exemple d'étude\n",
- "\n",
- "On s'intéresse (pour les exercices 1, 2 et 3) au problème de Cauchy\n",
- "\n",
- "$$\n",
- " \\dot{x}(t) = (1-2t) x(t), \\quad x(0) = x_0 = 1\n",
- "$$\n",
- "\n",
- "sur l'intervalle $[0, 3]$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Définition du problème de Cauchy\n",
- "f(t, x) = (1-2t)x # Second membre f(t, x)\n",
- "x0 = 1.0 # Condition initiale\n",
- "tspan = (0.0, 3.0); # Intervalle de temps"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Solution analytique\n",
- "function sol(t)\n",
- " return exp(t-t^2)\n",
- "end;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Estimation de la constante de Lipschitz de f sur [0, 3]\n",
- "# Voir Théorème 8.2.2 pour l'utilité de cette estimation\n",
- "function dfx(t)\n",
- " return 1-2t\n",
- "end \n",
- "L = maximum(abs.(dfx.(range(0, stop=3, length=1000))))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## La méthode d'Euler implicite\n",
- "\n",
- "La méthode d'Euler implicite est donnée par :\n",
- "\n",
- "$$\n",
- "\\left\\{\n",
- "\\begin{array}{l}\n",
- "x_{n+1} = x_n + h f(t_n + h, x_{n+1}) = G(x_{n+1}) \\\\\n",
- "x_0 = x(t_0)\n",
- "\\end{array}\n",
- "\\right.\n",
- "$$\n",
- "\n",
- "### Exercice 1\n",
- "\n",
- "1. Implémenter la méthode d'Euler implicite avec le point fixe (penser à voir le polycopié Section 8.2).\n",
- "2. Pourquoi si $h \\ge 0.2$, l'algorithme du point fixe peut s'arrêter au bout du nombre d'itérations max et donc ne pas converger pour la méthode d'Euler implicite ?\n",
- "3. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
- "4. Tracer l'erreur globale de la méthode d'Euler implicite en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$.\n",
- "\n",
- "**Attention** : pour l'algorithme du point fixe, faites attention aux critères d'arrêts (il y en a 2) ! Voir votre polycopié Section 8.2. Vous fixerez la valeur de la tolérance à $10^{-6}$ et le nombre maximum d'itérations à $1000$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "0.00447987015488227"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "# AJOUTER VOTRE CODE ICI\n",
- "\n",
- "# F(y) = y - G(y)\n",
- "\n",
- "G(y) = (y + sin(y))/3\n",
- "\n",
- "y = 1\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n",
- "y = G(y)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 11,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "0.0014932950464876392"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "F(y) = y - G(y)\n",
- "F(y)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## La méthode des trapèzes\n",
- "\n",
- "La méthode des trapèzes est donnée par le tableau de Butcher :\n",
- "\n",
- "$$\n",
- " \\begin{array}{c | c c}\n",
- " 0 & 0 & 0 \\\\[0.2em]\n",
- " 1 & 1/2 & 1/2 \\\\[0.2em]\n",
- " \\hline\n",
- " & 1/2 & 1/2 \\\\\n",
- " \\end{array}\n",
- "$$\n",
- "\n",
- "### Exercice 2\n",
- "\n",
- "1. Implémenter la méthode des trapèzes avec le point fixe.\n",
- "2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
- "3. Tracer l'erreur globale de la méthode des trapèzes. Quel est l'ordre de convergence de la méthode des trapèzes ?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# AJOUTER VOTRE CODE ICI"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## La méthode de Gauss à 2 étages\n",
- "\n",
- "La méthode de Gauss à 2 étages est donnée par le tableau de Butcher :\n",
- "\n",
- "$$\n",
- " \\begin{array}{c | c c}\n",
- " 1/2 - \\sqrt{3}/6 & 1/4 & 1/4 - \\sqrt{3}/6 \\\\[0.2em]\n",
- " 1/2 + \\sqrt{3}/6 & 1/4 + \\sqrt{3}/6 & 1/4 \\\\[0.2em]\n",
- " \\hline\n",
- " & 1/2 & 1/2 \\\\\n",
- " \\end{array}\n",
- "$$\n",
- "\n",
- "### Exercice 3\n",
- "\n",
- "1. Implémenter la méthode de Gauss à 2 étages avec le point fixe.\n",
- "2. Tracer la solution approchée et la solution exacte sur le même graphique pour différentes valeurs de $h$ que vous choisirez pour illustrer la convergence de la méthode.\n",
- "3. Tracer l'erreur globale de la méthode de Gauss à 2 étages. Quel est l'ordre de convergence de la méthode de Gauss à 2 étages ?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# AJOUTER VOTRE CODE ICI"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Un autre exemple\n",
- "\n",
- "On considère à partir de maintenant l'équation différentielle en dimension 2 :\n",
- "\n",
- "$$\n",
- " \\dot{x}_1(t) = x_2(t), \\quad \\dot{x}_2(t) = - x_1(t).\n",
- "$$\n",
- "\n",
- "On peut montrer facilement que la norme de $x(t) = (x_1(t), x_2(t))$ est constante le long des solutions :\n",
- "\n",
- "$$\n",
- " \\frac{\\mathrm{d}}{\\mathrm{d} t} \\|x(t)\\|^2 = 2\\, \\left( x(t) \\,|\\, \\dot{x}(t) \\right) = 2 \\left( x_1(t) x_2(t) - x_2(t) x_1(t) \\right) = 0.\n",
- "$$\n",
- "\n",
- "### Exercice 4\n",
- "\n",
- "On considère le problème de Cauchy associé de condition initiale $x_0 = (1, 0)$.\n",
- "\n",
- "1. Afficher l'approximation de la solution sur $[0, 10]$ pour les méthodes :\n",
- "- Euler explicite ;\n",
- "- Euler implicite ;\n",
- "- Trapèzes ;\n",
- "- Gauss à 2 étages.\n",
- "2. Commentaires.\n",
- "\n",
- "**Attention :** vous ferez un affichage dans le plan $(x_1, x_2)$. Vous fixerez le nombre de pas à $N=100$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# AJOUTER VOTRE CODE ICI"
- ]
- }
- ],
- "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/runge-kutta.ipynb b/tp/runge-kutta.ipynb
deleted file mode 100644
index 202f1ac46a337102128cd9f0144850d060b92288..0000000000000000000000000000000000000000
--- a/tp/runge-kutta.ipynb
+++ /dev/null
@@ -1,771 +0,0 @@
-{
- "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",
- "# Méthodes de Runge-Kutta\n",
- "\n",
- "- Date : 2023-2024\n",
- "- Durée approximative : 1h15"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "Introduction\n",
- "</span>\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "**Méthode à un pas explicite.** On désire calculer une approximation de la solution, notée $x(\\cdot, t_0, x_0)$, sur l'intervalle \n",
- "$I := [t_0, t_f]$ du problème de Cauchy\n",
- "\\begin{equation*}\n",
- " \\dot{x}(t) = f(t, x(t)), \\quad x(t_0) = x_0,\n",
- "\\end{equation*}\n",
- "où $f \\colon \\mathbb{R} \\times \\mathbb{R}^n \\to \\mathbb{R}^n$ est une application suffisamment régulière. \n",
- "On considère pour cela une subdivision de l'intervalle $I$ : \n",
- "$$\n",
- " t_0 < t_1 < \\cdots < t_N := t_f\n",
- "$$\n",
- "\n",
- "On note $h_i := t_{i+1} - t_i$ pour $i \\in \\llbracket 0, N-1 \\rrbracket$, les **pas** de la subdivision et \n",
- "$h_{\\max} := \\max_i(h_i)$ le pas le plus long. L'idée est de calculer une approximation de \n",
- "la solution en les points de discrétisation. On souhaite donc approcher $x(t_i, t_0, x_0)$ pour \n",
- "$i \\in \\llbracket 1, N \\rrbracket$. \n",
- "\n",
- "On appelle *méthode à un pas explicite*, toute méthode calculant un N-uplet $(x_1, \\ldots, x_N) \\in {(\\mathbb{R}^n)}^N$ tel que la valeur $x_{i+1}$ est calculée en fonction de $t_i$, $h_i$ et $x_i$ de la forme :\n",
- "\n",
- "$$\n",
- " x_{i+1} = x_i + h_i\\, \\Phi(t_i, x_i, h_i).\n",
- "$$\n",
- "\n",
- "Pour approcher la solution, nous utiliserons une méthode à un pas. Si on note $x_h^* := (x_1^*, \\ldots, x_N^*)$ une solution au problème discrétisé (obtenue par la méthode à un pas considérée), alors on souhaite que \n",
- "$$\n",
- " x_i^* \\approx x(t_i, t_0, x_0) \\quad \\text{pour tout } i \\in \\llbracket 1, N \\rrbracket.\n",
- "$$\n",
- "\n",
- "**Ordre de convergence.** On considère à partir de maintenant un pas de temps $h$ **uniforme**. Une méthode à un pas explicite est **convergente** si pour toute solution $x(\\cdot, t_0, x_0)$, la suite ${(x_i)}_i$ définie par $x_{i+1} = x_i + h\\, \\Phi(t_i, x_i, h)$ vérifie \n",
- "$$\n",
- " \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\to 0 \n",
- " \\quad\\text{quand}\\quad h \\to 0.\n",
- "$$\n",
- "\n",
- "Si la convergence est d'ordre $p$ alors il existe une constante $C \\ge 0$ telle que\n",
- "\n",
- "$$\n",
- " E := \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\le C\\, h^p.\n",
- "$$\n",
- "\n",
- "Faisons l'hypothèse que $E = M\\, h^p$ pour un certain $M \\ge 0$. En passant au logarithme, on obtient\n",
- "\n",
- "$$\n",
- " \\log(E) = \\log(M) + p\\, \\log(h).\n",
- "$$\n",
- "\n",
- "Nous en déduisons que si on trace $\\log(E)$ en fonction de $\\log(h)$, on doit obtenir une droite de pente $p$. C'est ce que nous allons vérifier dans ce TP.\n",
- "\n",
- "**Objectif.** Nous allons dans ce TP, implémenter quelques méthodes (à un pas) de [Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta) explicites et \n",
- "étudier leur convergence. On appelle *méthode de Runge-Kutta explicite à $s$ étages*, la méthode définie par le schéma\n",
- "\n",
- "$$\n",
- " \\begin{aligned}\n",
- " k_1 &= f(t_0, x_0) \\\\\n",
- " k_2 &= f(t_0 + c_2 h, x_0 + h a_{21} k_1) \\\\\n",
- " \\vdots & \\\\[-1em]\n",
- " k_s &= f(t_0 + c_s h, x_0 + h \\sum_{j=1}^{s-1} a_{sj} k_j) \\\\\n",
- " x_1 &= x_0 + h \\sum_{i=1}^{s} b_i k_i,\n",
- " \\end{aligned}\n",
- "$$\n",
- "où les coefficients $c_i$, $a_{ij}$ et $b_i$ sont des constantes qui définissent précisément le schéma. \n",
- "On supposera toujours dans la suite que $c_1 = 0$ et $c_i = \\sum_{j=1}^{i-1} a_{ij}$, pour $i = 2, \\ldots, s$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# load packages\n",
- "using Plots\n",
- "using Plots.PlotMeasures\n",
- "using Polynomials\n",
- "\n",
- "#\n",
- "px = PlotMeasures.px;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Fonctons auxiliaires \n",
- "\n",
- "function method_infos(method)\n",
- "\n",
- " if method == :euler \n",
- " method_func = euler\n",
- " method_name = \"Euler\"\n",
- " methode_stages = 1\n",
- " elseif method == :runge \n",
- " method_func = runge\n",
- " method_name = \"Runge\"\n",
- " methode_stages = 2\n",
- " elseif method == :heun\n",
- " method_func = heun\n",
- " method_name = \"Heun\"\n",
- " methode_stages = 3\n",
- " elseif method == :rk4\n",
- " method_func = rk4\n",
- " method_name = \"RK4\"\n",
- " methode_stages = 4\n",
- " else \n",
- " error(\"Méthode d'intégration non reconnue\")\n",
- " end\n",
- "\n",
- " return method_func, method_name, methode_stages\n",
- "\n",
- "end\n",
- "\n",
- "function convergence(method, f, x0, tspan, sol, Nspan)\n",
- "\n",
- " # Récupération des informations sur la méthode\n",
- " method_func, method_name, methode_stages = method_infos(method)\n",
- " \n",
- " # Ecriture des choix de paramètres\n",
- " println(\"Méthode d'intégration : \", method_name)\n",
- "\n",
- " plts = [] # Liste des graphiques\n",
- "\n",
- " for N ∈ Nspan\n",
- "\n",
- " # Solution approchée\n",
- " ts, xs = method_func(f, x0, tspan, N) # Appel de la méthode d'intégraton\n",
- "\n",
- " # Affichage de la solution approchée\n",
- " plt = plot(ts, xs, label=method_name, marker=:circle)\n",
- "\n",
- " # Affichage de la solution analytique\n",
- " plot!(plt, sol, label=\"Analytic\")\n",
- "\n",
- " # Mise en forme du graphique\n",
- " plot!(plt, xlims=(tspan[1]-0.1, tspan[2]+0.1), ylims=(-0.2, 1.65), \n",
- " title=\"h=$(round((tspan[2]-tspan[1])/N, digits=2))\", titlefont = font(12,\"Calibri\"),\n",
- " xlabel=\"t\", ylabel=\"x(t)\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
- "\n",
- " # Ajout du graphique à la liste des graphiques\n",
- " push!(plts, plt)\n",
- "\n",
- " end\n",
- "\n",
- " return plts\n",
- "\n",
- "end\n",
- "\n",
- "xlims_ = (1e-4, 1e0)\n",
- "ylims_ = (1e-13, 1e1)\n",
- "xlims_nfe_ = (1e1, 1e4)\n",
- "ylims_nfe_ = (1e-13, 1e1)\n",
- "\n",
- "function ordre(method, f, x0, tspan, sol, hspan, nfespan; \n",
- " xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)\n",
- "\n",
- " plts = [] # Liste des graphiques\n",
- "\n",
- " #\n",
- " plt1 = plot(xaxis=:log, yaxis=:log); push!(plts, plt1)\n",
- " plt2 = plot(xaxis=:log, yaxis=:log); push!(plts, plt2)\n",
- "\n",
- " #\n",
- " plts = ordre(plts, method, f, x0, tspan, sol, hspan, nfespan, \n",
- " xlims=xlims, ylims=ylims, xlims_nfe=xlims_nfe, ylims_nfe=ylims_nfe)\n",
- "\n",
- " # Mise en forme des graphiques\n",
- " plot!(plts[1], titlefont = font(12, \"Calibri\"), legend=:topleft,\n",
- " xlabel=\"h\", ylabel=\"Error\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
- "\n",
- " plot!(plts[2], titlefont = font(12, \"Calibri\"), legend=:topright,\n",
- " xlabel=\"Calls to f\", ylabel=\"Error\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
- "\n",
- " return plts\n",
- " \n",
- "end\n",
- "\n",
- "function ordre(plts_in, method, f, x0, tspan, sol, hspan, nfespan; \n",
- " xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)\n",
- "\n",
- " # Copie du graphique d'entrée\n",
- " plts = deepcopy(plts_in)\n",
- "\n",
- " # Récupération des informations sur la méthode\n",
- " method_func, method_name, methode_stages = method_infos(method)\n",
- " \n",
- " # Ecriture des choix de paramètres\n",
- " println(\"Méthode d'intégration : \", method_name)\n",
- "\n",
- " # Les différents nombre de pas de temps\n",
- " Nspan = round.(Int, (tspan[2]-tspan[1]) ./ hspan)\n",
- "\n",
- " # Calcul de l'erreur\n",
- " err = []\n",
- "\n",
- " for N ∈ Nspan\n",
- "\n",
- " # Solution approchée\n",
- " ts, xs = method_func(f, x0, tspan, N)\n",
- "\n",
- " # On calcule l'erreur en norme infinie\n",
- " push!(err, maximum(abs.(xs .- sol.(ts))))\n",
- " \n",
- " end\n",
- "\n",
- " # calcul par régression linéaire de la pente de la droite et de l'ordonnée à l'origine\n",
- " reg = fit(log10.(hspan), log10.(err), 1)\n",
- " K = 10^reg[0]\n",
- " p = reg[1]\n",
- " println(\"\\nconstante du grand O : K = $(round(K, digits=5))\")\n",
- " println(\"ordre de convergence : p = $(round(p, digits=5))\")\n",
- "\n",
- " # Affichage de l'erreur en fonction du pas de temps: on enlève la constante K\n",
- " plot!(plts[1], hspan, err, xaxis=:log, yaxis=:log, label=\"$method_name\", marker=:circle)\n",
- "\n",
- " # Affichage de la droite de régression\n",
- " #plot!(plt, hspan, hspan .^ p, label=\"$method_name Regression\", linestyle=:dash)\n",
- " \n",
- " # Mise en forme du graphique\n",
- " plot!(plts[1], xlims=xlims, ylims=ylims)\n",
- "\n",
- " # Calcul de l'erreur en fonction du nombre d'appels à la fonction f\n",
- " err = []\n",
- "\n",
- " for Nfe ∈ Nfespan\n",
- "\n",
- " #\n",
- " N = round(Int, Nfe / methode_stages)\n",
- "\n",
- " # Solution approchée\n",
- " ts, xs = method_func(f, x0, tspan, N)\n",
- "\n",
- " # On calcule l'erreur en norme infinie\n",
- " push!(err, maximum(abs.(xs .- sol.(ts))))\n",
- " \n",
- " end\n",
- "\n",
- " # Affchage de l'erreur en fonction du nombre d'appels à la fonction f\n",
- " # Nfespan = methode_stages .* Nspan\n",
- " plot!(plts[2], Nfespan, err, xaxis=:log, yaxis=:log, label=\"$method_name\", marker=:circle)\n",
- "\n",
- " # Mise en forme du graphique\n",
- " plot!(plts[2], xlims=xlims_nfe, ylims=ylims_nfe)\n",
- "\n",
- " return plts\n",
- " \n",
- "end;\n",
- "\n",
- "hspan_ = 10 .^ range(-4, stop=0, length=20)\n",
- "Nfespan_ = 10 .^ range( 0, stop=4, length=20);"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "L'exemple d'étude\n",
- "</span>"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "On s'intéresse dans ce TP au problème de Cauchy\n",
- "\n",
- "$$\n",
- " \\dot{x}(t) = (1-2t) x(t), \\quad x(0) = x_0 = 1\n",
- "$$\n",
- "\n",
- "sur l'intervalle $[0, 3]$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Définition du problème de Cauchy\n",
- "f(t, x) = (1-2t)x # Second membre f(t, x)\n",
- "x0 = 1.0 # Condition initiale\n",
- "tspan = (0.0, 3.0); # Intervalle de temps"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Exercice 1\n",
- "\n",
- "* Calculer la solution analytique de l'exemple ci-dessus et la coder dans la fonction `sol` ci-dessous."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Solution analytique\n",
- "function sol(t)\n",
- " return 0 # TO UPDATE\n",
- "end;"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "La méthode d'Euler\n",
- "</span>"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "La méthode d'Euler est donnée par :\n",
- "\n",
- "$$\n",
- "\\left\\{\n",
- "\\begin{array}{l}\n",
- "x_{n+1} = x_n + h f(t_n, x_n) \\\\\n",
- "x_0 = x(t_0)\n",
- "\\end{array}\n",
- "\\right.\n",
- "$$\n",
- "\n",
- "### Exercice 2\n",
- "\n",
- "1. Implémenter la méthode d'Euler dans la fonction `euler` ci-dessous.\n",
- "2. Vérifier la convergence de la méthode d'Euler sur l'exemple donné dans la cellule d'après."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Implémentation de la méthode d'Euler\n",
- "function euler(f, x0, tspan, N)\n",
- " t0, tf = tspan\n",
- " h = (tf - t0) / N\n",
- " t = t0\n",
- " x = x0\n",
- " ts = [t0]\n",
- " xs = [x0]\n",
- " for i in 1:N\n",
- " x = x # TO UPDATE\n",
- " t = t # TO UPDATE\n",
- " push!(ts, t)\n",
- " push!(xs, x)\n",
- " end\n",
- " return ts, xs\n",
- "end;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Convergence de la méthode sur l'exemple\n",
- "method = :euler\n",
- "Nspan = [5, 10, 20]\n",
- "\n",
- "# Calcul des graphiques\n",
- "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
- "\n",
- "# Affichage des graphiques\n",
- "plot(plts..., layout=(1, length(Nspan)), size=(800, 400))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "3. Afficher l'erreur globale de la méthode d'Euler en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Ordre de convergence de la méthode\n",
- "method = :euler\n",
- "hspan = hspan_[1e-7 .≤ hspan_ .≤ 1]\n",
- "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e7]\n",
- "\n",
- "# Calcul du graphique\n",
- "plt_order_euler = ordre(method, f, x0, tspan, sol, hspan, Nfespan)\n",
- "\n",
- "# Affichage du graphique\n",
- "plot(plt_order_euler..., layout=(1, 2), size=(800, 400))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "La méthode de Runge\n",
- "</span>"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "La méthode de Runge est donnée par :\n",
- "\n",
- "$$\n",
- "\\left\\{\n",
- "\\begin{array}{l}\n",
- "\\displaystyle x_{n+1} = x_n + h\\, f\\left(t_{n} + \\frac{h}{2}, ~x_n + \\frac{h}{2} f(t_n, x_n)\\right) \\\\\n",
- "x_0 = x(t_0)\n",
- "\\end{array}\n",
- "\\right.\n",
- "$$\n",
- "\n",
- "### Exercice 3\n",
- "\n",
- "1. Implémenter la méthode de Runge dans la fonction `runge` ci-dessous.\n",
- "2. Vérifier la convergence de la méthode de Runge sur l'exemple donné dans la cellule d'après.\n",
- "3. Calculer l'erreur globale de la méthode de Runge en fonction de $h$ et vérifier que l'erreur est bien en $O(h^2)$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function runge(f, x0, tspan, N)\n",
- " t0, tf = tspan\n",
- " h = (tf - t0) / N\n",
- " t = t0\n",
- " x = x0\n",
- " ts = [t0]\n",
- " xs = [x0]\n",
- " for i in 1:N\n",
- " k1 = f(t, x)\n",
- " k2 = 0 # TO UPDATE\n",
- " x = x # TO UPDATE\n",
- " t = t # TO UPDATE\n",
- " push!(ts, t)\n",
- " push!(xs, x)\n",
- " end\n",
- " return ts, xs\n",
- "end;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Convergence de la méthode sur l'exemple\n",
- "method = :runge\n",
- "Nspan = [5, 10, 20]\n",
- "\n",
- "# Calcul des graphiques\n",
- "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
- "\n",
- "# Affichage des graphiques\n",
- "plot(plts..., layout=(1, length(Nspan)), size=(800, 400))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Ordre de convergence de la méthode\n",
- "method = :runge\n",
- "hspan = hspan_[1e-5 .≤ hspan_ .≤ 1]\n",
- "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e6]\n",
- "\n",
- "# Calcul du graphique\n",
- "plt_order_runge = ordre(plt_order_euler, method, f, x0, tspan, sol, hspan, Nfespan)\n",
- "\n",
- "# Affichage du graphique\n",
- "plot(plt_order_runge..., layout=(1, 2), size=(800, 400))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "La méthode de Heun\n",
- "</span>"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "La méthode de Heun est donnée par :\n",
- "\n",
- "$$\n",
- "\\left\\{\n",
- "\\begin{array}{l}\n",
- " \\displaystyle k_1 = f(t_n, x_n) \\\\[1em]\n",
- " \\displaystyle k_2 = f(t_n + \\frac{h}{3}, x_n + \\frac{h}{3} k_1) \\\\[1em]\n",
- " \\displaystyle k_3 = f(t_n + \\frac{2}{3}h, x_n + \\frac{2}{3}h\\, k_2) \\\\[1em]\n",
- " \\displaystyle x_{n+1} = x_n + h\\left(\\frac{1}{4} k_1 + \\frac{3}{4} k_3\\right), \\quad \n",
- " \\displaystyle x_0 = x(t_0)\n",
- "\\end{array}\n",
- "\\right.\n",
- "$$\n",
- "\n",
- "### Exercice 4\n",
- "\n",
- "1. Implémenter la méthode de Heun dans la fonction `heun` ci-dessous.\n",
- "2. Vérifier la convergence de la méthode de Heun sur l'exemple donné dans la cellule d'après.\n",
- "3. Calculer l'erreur globale de la méthode de Heun en fonction de $h$ et vérifier que l'erreur est bien en $O(h^3)$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function heun(f, x0, tspan, N)\n",
- " t0, tf = tspan\n",
- " h = (tf - t0) / N\n",
- " t = t0\n",
- " x = x0\n",
- " ts = [t0]\n",
- " xs = [x0]\n",
- " for i in 1:N\n",
- " k1 = f(t, x)\n",
- " k2 = 0 # TO UPDATE\n",
- " k3 = 0 # TO UPDATE \n",
- " x = x # TO UPDATE\n",
- " t = t # TO UPDATE\n",
- " push!(ts, t)\n",
- " push!(xs, x)\n",
- " end\n",
- " return ts, xs\n",
- "end;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Convergence de la méthode sur l'exemple\n",
- "method = :heun\n",
- "Nspan = [5, 10, 20]\n",
- "\n",
- "# Calcul des graphiques\n",
- "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
- "\n",
- "# Affichage des graphiques\n",
- "plot(plts..., layout=(1, length(Nspan)), size=(800, 400))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Ordre de convergence de la méthode\n",
- "method = :heun\n",
- "hspan = hspan_[1e-4 .≤ hspan_ .≤ 1]\n",
- "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e5]\n",
- "\n",
- "# Calcul du graphique\n",
- "plt_order_heun = ordre(plt_order_runge, method, f, x0, tspan, sol, hspan, Nfespan)\n",
- "\n",
- "# Affichage du graphique\n",
- "plot(plt_order_heun..., layout=(1, 2), size=(800, 400))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "<span style=\"margin:0px;\n",
- " padding:8px;\n",
- " background-color:#afa;\n",
- " border:2px solid #bbffbb;\n",
- " border-radius:20px;\n",
- " font-weight:bold;\n",
- " font-size:1.5em;\n",
- " text-align:center;\">\n",
- "La méthode de Runge-Kutta d'ordre 4\n",
- "</span>"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "La méthode de Runge-Kutta d'ordre 4 est donnée par :\n",
- "\n",
- "$$\n",
- "\\left\\{\n",
- "\\begin{array}{l}\n",
- " \\displaystyle k_1 = f(t_n, x_n) \\\\[1em]\n",
- " \\displaystyle k_2 = f(t_n + \\frac{h}{2}, x_n + \\frac{h}{2} k_1) \\\\[1em]\n",
- " \\displaystyle k_3 = f(t_n + \\frac{h}{2}, x_n + \\frac{h}{2} k_2) \\\\[1em]\n",
- " \\displaystyle k_4 = f(t_n + h, x_n + h\\, k_3) \\\\[1em]\n",
- " \\displaystyle x_{n+1} = x_n + \\frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4), \\quad \n",
- " \\displaystyle x_0 = x(t_0)\n",
- "\\end{array}\n",
- "\\right.\n",
- "$$\n",
- "\n",
- "### Exercice 5\n",
- "\n",
- "1. Implémenter la méthode de Runge-Kutta d'ordre 4 dans la fonction `rk4` ci-dessous.\n",
- "2. Vérifier la convergence de la méthode de Runge-Kutta d'ordre 4 sur l'exemple donné dans la cellule d'après.\n",
- "3. Calculer l'erreur globale de la méthode de Runge-Kutta d'ordre 4 en fonction de $h$ et vérifier que l'erreur est bien en $O(h^4)$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "function rk4(f, x0, tspan, N)\n",
- " t0, tf = tspan\n",
- " h = (tf - t0) / N\n",
- " t = t0\n",
- " x = x0\n",
- " ts = [t0]\n",
- " xs = [x0]\n",
- " for i in 1:N\n",
- " k1 = f(t, x)\n",
- " k2 = 0 # TO UPDATE\n",
- " k3 = 0 # TO UPDATE\n",
- " k4 = 0 # TO UPDATE\n",
- " x = x # TO UPDATE\n",
- " t = t # TO UPDATE\n",
- " push!(ts, t)\n",
- " push!(xs, x)\n",
- " end\n",
- " return ts, xs\n",
- "end;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Convergence de la méthode sur l'exemple\n",
- "method = :rk4\n",
- "Nspan = [5, 10, 20]\n",
- "\n",
- "# Calcul des graphiques\n",
- "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
- "\n",
- "# Affichage des graphiques\n",
- "plot(plts..., layout=(1, length(Nspan)), size=(800, 400))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Ordre de convergence de la méthode\n",
- "method = :rk4\n",
- "hspan = hspan_[1e-3 .≤ hspan_ .≤ 1]\n",
- "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e4]\n",
- "\n",
- "# Calcul du graphique\n",
- "plt_order_rk4 = ordre(plt_order_heun, method, f, x0, tspan, sol, hspan, Nfespan)\n",
- "\n",
- "# Affichage du graphique\n",
- "plot(plt_order_rk4..., layout=(1, 2), size=(800, 400))"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.10.3",
- "language": "julia",
- "name": "julia-1.10"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.10.3"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/tp/simple-shooting.ipynb b/tp/simple-shooting.ipynb
deleted file mode 100644
index 95eb36385c58c711dcad707cc724d88a5088bdbc..0000000000000000000000000000000000000000
--- a/tp/simple-shooting.ipynb
+++ /dev/null
@@ -1,368 +0,0 @@
-{
- "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
-}
diff --git a/tp/space.jl b/tp/space.jl
deleted file mode 100644
index 1f9c72aaf5b65ac8ff8f09316865a8daefccd9ee..0000000000000000000000000000000000000000
--- a/tp/space.jl
+++ /dev/null
@@ -1,849 +0,0 @@
-module Space
-
- using Plots
- using Plots.PlotMeasures
- using SparseArrays
- using Printf
- using DifferentialEquations
- using LinearAlgebra # norm
-
- export animation
-
- # --------------------------------------------------------------------------------------------
- x0 = [-42272.67, 0, 0, -5796.72] # état initial
- μ = 5.1658620912*1e12
- rf = 42165.0
- t0 = 0
- global γ_max = nothing
- global F_max = nothing
-
- # Maximizing control
- function control(p)
- u = zeros(eltype(p),2)
- u[1] = p[3]*γ_max/norm(p[3:4])
- u[2] = p[4]*γ_max/norm(p[3:4])
- return u
- end
-
- # --------------------------------------------------------------------------------------------
- # Données pour la trajectoire durant le transfert
- mutable struct Transfert
- initial_adjoint
- duration
- flow
- end
-
- # données pour la trajectoire pré-transfert
- mutable struct PreTransfert
- initial_point
- duration
- end
-
- # données pour la trajectoire post-transfert
- mutable struct PostTransfert
- duration
- end
-
- # --------------------------------------------------------------------------------------------
- # Default options for flows
- # --------------------------------------------------------------------------------------------
- function __abstol()
- return 1e-10
- end
-
- function __reltol()
- return 1e-10
- end
-
- function __saveat()
- return []
- end
-
- # --------------------------------------------------------------------------------------------
- #
- # Vector Field
- #
- # --------------------------------------------------------------------------------------------
- struct VectorField f::Function end
-
- function (vf::VectorField)(x) # https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
- return vf.f(x)
- end
-
- # Flow of a vector field
- function Flow(vf::VectorField)
-
- function rhs!(dx, x, dummy, t)
- dx[:] = vf(x)
- end
-
- function f(tspan, x0; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
- ode = ODEProblem(rhs!, x0, tspan)
- sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
- return sol
- end
-
- function f(t0, x0, t; abstol=__abstol(), reltol=__reltol(), saveat=__saveat())
- sol = f((t0, t), x0; abstol=abstol, reltol=reltol, saveat=saveat)
- n = size(x0, 1)
- return sol[1:n, end]
- end
-
- return f
-
- end
-
- # Kepler's law
- function kepler(x)
- n = size(x, 1)
- dx = zeros(eltype(x), n)
- x1 = x[1]
- x2 = x[2]
- x3 = x[3]
- x4 = x[4]
- dsat = norm(x[1:2]); r3 = dsat^3;
- dx[1] = x3
- dx[2] = x4
- dx[3] = -μ*x1/r3
- dx[4] = -μ*x2/r3
- return dx
- end
-
- KeplerFlow = Flow(VectorField(kepler));
-
- function animation(p0, tf, f, γ_max; fps=10, nFrame=100)
- #
- t_pre_transfert = 2*5.302395297743802
- x_pre_transfert = x0
- pre_transfert_data = PreTransfert(x_pre_transfert, t_pre_transfert)
- #
- transfert_data = Transfert(p0, tf, f)
- #
- t_post_transfert = 20
- post_transfert_data = PostTransfert(t_post_transfert)
- #
- __animation(pre_transfert_data, transfert_data, post_transfert_data, KeplerFlow, γ_max, fps=fps, nFrame=nFrame)
- end
-
- function __animation(pre_transfert_data, transfert_data, post_transfert_data, f0, _γ_max; fps=10, nFrame=200)
-
- #
- global F_max = _γ_max*2000*10^3/3600^2
- global γ_max = _γ_max
-
- # pré-transfert
- t_pre_transfert = pre_transfert_data.duration
- x_pre_transfert = pre_transfert_data.initial_point
-
- # transfert
- p0_transfert = transfert_data.initial_adjoint
- tf_transfert = transfert_data.duration
- f = transfert_data.flow
-
- # post-trasfert
- t_post_transfert = post_transfert_data.duration
-
- # On calcule les orbites initiale et finale
- r0 = norm(x0[1:2])
- v0 = norm(x0[3:4])
- a = 1.0/(2.0/r0-v0*v0/μ)
- t1 = r0*v0*v0/μ - 1.0;
- t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ);
- e_ellipse = norm([t1 t2])
- p_orb = a*(1-e_ellipse^2);
- n_theta = 151
- Theta = range(0.0, stop=2*pi, length=n_theta)
- X1_orb_init = zeros(n_theta)
- X2_orb_init = zeros(n_theta)
- X1_orb_arr = zeros(n_theta)
- X2_orb_arr = zeros(n_theta)
-
- for i in 1:n_theta
- theta = Theta[i]
- r_orb = p_orb/(1+e_ellipse*cos(theta));
- # Orbite initiale
- X1_orb_init[i] = r_orb*cos(theta);
- X2_orb_init[i] = r_orb*sin(theta);
- # Orbite finale
- X1_orb_arr[i] = rf*cos(theta) ;
- X2_orb_arr[i] = rf*sin(theta);
- end
-
- # Centre de la fenêtre
- cx = 40000
- cy = -7000
-
- # Taille de la fenêtre
- w = 1600
- h = 900
- ρ = h/w
-
- # Limites de la fenêtre
- ee = 0.5
- xmin = minimum([X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);
- xmax = maximum([X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);
- ymin = minimum([X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);
- ymax = maximum([X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);
-
- Δy = ymax-ymin
- Δx = Δy/ρ
- Δn = Δx - (xmax-xmin)
- xmin = xmin - Δn/2
- xmax = xmax + Δn/2
-
- # Trajectoire pré-transfert
- traj_pre_transfert = f0((0.0, t_pre_transfert), x_pre_transfert)
- times_pre_transfert = traj_pre_transfert.t
- n_pre_transfert = size(times_pre_transfert, 1)
- x1_pre_transfert = [traj_pre_transfert[1, j] for j in 1:n_pre_transfert ]
- x2_pre_transfert = [traj_pre_transfert[2, j] for j in 1:n_pre_transfert ]
- v1_pre_transfert = [traj_pre_transfert[3, j] for j in 1:n_pre_transfert ]
- v2_pre_transfert = [traj_pre_transfert[4, j] for j in 1:n_pre_transfert ]
-
- # Trajectoire pendant le transfert
- traj_transfert = f((t0, tf_transfert), x0, p0_transfert)
- times_transfert = traj_transfert.t
- n_transfert = size(times_transfert, 1)
- x1_transfert = [traj_transfert[1, j] for j in 1:n_transfert ]
- x2_transfert = [traj_transfert[2, j] for j in 1:n_transfert ]
- v1_transfert = [traj_transfert[3, j] for j in 1:n_transfert ]
- v2_transfert = [traj_transfert[4, j] for j in 1:n_transfert ]
- u_transfert = zeros(2, length(times_transfert))
- for j in 1:n_transfert
- u_transfert[:,j] = control(traj_transfert[5:8, j])
- end
- u_transfert = u_transfert./γ_max # ||u|| ≤ 1
-
- # post-transfert
- x_post_transfert = traj_transfert[:,end]
-
- # Trajectoire post-transfert
- traj_post_transfert = f0((0.0, t_post_transfert), x_post_transfert)
- times_post_transfert = traj_post_transfert.t
- n_post_transfert = size(times_post_transfert, 1)
- x1_post_transfert = [traj_post_transfert[1, j] for j in 1:n_post_transfert ]
- x2_post_transfert = [traj_post_transfert[2, j] for j in 1:n_post_transfert ]
- v1_post_transfert = [traj_post_transfert[3, j] for j in 1:n_post_transfert ]
- v2_post_transfert = [traj_post_transfert[4, j] for j in 1:n_post_transfert ]
-
- # Angles de rotation du satellite pendant le pré-transfert
- # Et poussée normalisée entre 0 et 1
- θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])
- θ_pre_transfert = range(π/2, mod(θ0, 2*π), length=n_pre_transfert)
- F_pre_transfert = zeros(1, n_pre_transfert)
-
- # Angles de rotation du satellite pendant le transfert
- θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])
- F_transfert = zeros(1, n_transfert)
- for j in 1:n_transfert
- F_transfert[j] = norm(u_transfert[:,j])
- end
-
- # Angles de rotation du satellite pendant le post-transfert
- θ1 = atan(u_transfert[2, end], u_transfert[1, end])
- θ2 = atan(-x2_post_transfert[end], -x1_post_transfert[end])
- θ_post_transfert = range(mod(θ1, 2*π), mod(θ2, 2*π), length=n_post_transfert)
- F_post_transfert = zeros(1, n_post_transfert)
-
- # Etoiles
- Δx = xmax-xmin
- Δy = ymax-ymin
- ρ = Δy/Δx
- S = stars(ρ)
-
- # nombre total d'images
- nFrame = min(nFrame, n_pre_transfert+n_transfert+n_post_transfert);
-
- # Pour l'affichage de la trajectoire globale
- times = [times_pre_transfert[1:end-1];
- times_pre_transfert[end].+times_transfert[1:end-1];
- times_pre_transfert[end].+times_transfert[end].+times_post_transfert[1:end]]
- x1 = [x1_pre_transfert[1:end-1]; x1_transfert[1:end-1]; x1_post_transfert[:]]
- x2 = [x2_pre_transfert[1:end-1]; x2_transfert[1:end-1]; x2_post_transfert[:]]
- v1 = [v1_pre_transfert[1:end-1]; v1_transfert[1:end-1]; v1_post_transfert[:]]
- v2 = [v2_pre_transfert[1:end-1]; v2_transfert[1:end-1]; v2_post_transfert[:]]
- θ = [ θ_pre_transfert[1:end-1]; θ_transfert[1:end-1]; θ_post_transfert[:]]
- F = [ F_pre_transfert[1:end-1]; F_transfert[1:end-1]; F_post_transfert[:]]
-
- # plot thrust on/off
- th = [BitArray(zeros(n_pre_transfert-1));
- BitArray(ones(n_transfert-1));
- BitArray(zeros(n_post_transfert))]
-
- # plot trajectory
- pt = [BitArray(ones(n_pre_transfert-1));
- BitArray(ones(n_transfert-1));
- BitArray(zeros(n_post_transfert))]
-
- # Contrôle sur la trajectoire totale
- u_total = hcat([zeros(2, n_pre_transfert-1),
- u_transfert[:, 1:n_transfert-1],
- zeros(2, n_post_transfert)]...)
-
- # temps total
- temps_transfert_global = times[end]
-
- # pas de temps pour le transfert global
- if nFrame>1
- Δt = temps_transfert_global/(nFrame-1)
- else
- Δt = 0.0
- end
-
- # opacités des orbites initiale et finale
- op_initi = [range(0.0, 1.0, length=n_pre_transfert-1);
- range(1.0, 0.0, length=n_transfert-1);
- zeros(n_post_transfert)]
- op_final = [zeros(n_pre_transfert-1);
- range(0.0, 1.0, length=n_transfert-1);
- range(1.0, 0.0, length=int(n_post_transfert/4));
- zeros(n_post_transfert-int(n_post_transfert/4))]
-
- println("Fps = ", fps)
- println("NFrame = ", nFrame)
- println("Vitesse = ", nFrame/(temps_transfert_global*fps))
- println("Durée totale de la mission = ", temps_transfert_global, " h")
-
- # animation
- anim = @animate for i ∈ 1:nFrame
-
- # Δt : pas de temps
- # time_current : temps courant de la mission totale à l'itération i
- # i_current : indice tel que times[i_current] = time_current
- # w, h : width, height de la fenêtre
- # xmin, xmax, ymin, ymax : limites des axes du plot principal
- # X1_orb_init, X2_orb_init : coordonnées de l'orbite initiale
- # X1_orb_arr, X2_orb_arr : coordonnées de l'orbite finale
- # cx, cy : coordonées du centre de l'affichage du tranfert
- # S : data pour les étoiles
- # Δx, Δy : xmax-xmin, ymax-ymin
- # times : tous les temps de la mission complète, ie pre-transfert, transfert et post-transfert
- # x1, x2 : vecteur de positions du satellite
- # θ : vecteur d'orientations du satellite
- # th : vecteur de booléens - thrust on/off
- # u_total : vecteur de contrôles pour toute la mission
- # F_max, γ_max : poussée max
- # subplot_current : valeur du subplot courant
- # cam_x, cam_y : position de la caméra
- # cam_zoom : zoom de la caméra
-
- cam_x = cx
- cam_y = cy
- cam_zoom = 1
-
- time_current = (i-1)*Δt
- i_current = argmin(abs.(times.-time_current))
-
- px = background(w, h, xmin, xmax, ymin, ymax,
- X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
- cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom,
- op_initi[i_current], op_final[i_current], times, time_current)
-
- trajectoire!(px, times, x1, x2, θ, F, th, time_current, cx, cy, pt)
-
- subplot_current = 2
- subplot_current = panneau_control!(px, time_current, times, u_total,
- F_max, subplot_current)
-
- subplot_current = panneau_information!(px, subplot_current, time_current,
- i_current, x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init,
- X2_orb_init, X1_orb_arr, X2_orb_arr)
-
- end
-
- # enregistrement
- #gif(anim, "transfert-temps-min-original.mp4", fps=fps);
- gif(anim, "transfert-temps-min.gif", fps=fps);
-
- end;
-
- # --------------------------------------------------------------------------------------------
- # Preliminaries
- int(x) = floor(Int, x);
- rayon_Terre = 6371;
- sc_sat = 1000; # scaling for the satellite
-
- # --------------------------------------------------------------------------------------------
- #
- # Satellite
- #
- # --------------------------------------------------------------------------------------------
- # Corps du satellite
- @userplot SatBody
- @recipe function f(cp::SatBody)
- x, y = cp.args
- seriestype --> :shape
- fillcolor --> :goldenrod1
- linecolor --> :goldenrod1
- x, y
- end
-
- @userplot SatLine
- @recipe function f(cp::SatLine)
- x, y = cp.args
- linecolor --> :black
- x, y
- end
-
- # Bras du satellite
- @userplot SatBras
- @recipe function f(cp::SatBras)
- x, y = cp.args
- seriestype --> :shape
- fillcolor --> :goldenrod1
- linecolor --> :white
- x, y
- end
-
- # panneau
- @userplot PanBody
- @recipe function f(cp::PanBody)
- x, y = cp.args
- seriestype --> :shape
- fillcolor --> :dodgerblue4
- linecolor --> :black
- x, y
- end
-
- # flamme
- @userplot Flamme
- @recipe function f(cp::Flamme)
- x, y = cp.args
- seriestype --> :shape
- fillcolor --> :darkorange1
- linecolor --> false
- x, y
- end
-
- function satellite!(pl; subplot=1, thrust=false, position=[0;0], scale=1, rotate=0, thrust_val=1.0)
-
- # Fonctions utiles
- R(θ) = [ cos(θ) -sin(θ)
- sin(θ) cos(θ)] # rotation
- T(x, v) = x.+v # translation
- H(λ, x) = λ.*x # homotéthie
- SA(x, θ, c) = c .+ 2*[cos(θ);sin(θ)]'*(x.-c).*[cos(θ);sin(θ)]-(x.-c) # symétrie axiale
-
- #
- O = position
- Rθ = R(rotate)
-
- # Paramètres
- α = π/10.0 # angle du tube, par rapport à l'horizontal
- β = π/2-π/40 # angle des bras, par rapport à l'horizontal
- Lp = scale.*1.4*cos(α) # longueur d'un panneau
- lp = scale.*2.6*sin(α) # largueur d'un panneau
-
- # Param bras du satellite
- lb = scale.*3*cos(α) # longueur des bras
- eb = scale.*cos(α)/30 # demi largeur des bras
- xb = 0.0
- yb = scale.*sin(α)
-
- # Paramètres corps du satellite
- t = range(-α, α, length=50)
- x = scale.*cos.(t);
- y = scale.*sin.(t);
- Δx = scale.*cos(α)
- Δy = scale.*sin(α)
-
- # Paramètres flamme
- hF = yb # petite hauteur
- HF = 2*hF # grande hauteur
- LF = 5.5*Δx # longueur
-
- # Dessin bras du satellite
- M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
- [yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb]'], O);
- satbras!(pl, M[1,:], M[2,:], subplot=subplot)
- M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
- -[yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb']'], O);
- satbras!(pl, M[1,:], M[2,:], subplot=subplot)
-
- # Dessin corps du satellite
- M = T(Rθ*[ x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
- M = T(Rθ*[-x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
- M = T(Rθ*[scale.*[-cos(α), cos(α), cos(α), -cos(α)]'
- scale.*[-sin(α), -sin(α), sin(α), sin(α)]'], O);
- satbody!(pl, M[1,:], M[2,:], subplot=subplot) # interieur
-
- M = T(Rθ*[ x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
- M = T(Rθ*[-x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
- M = T(Rθ*[ x'.-2*Δx; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche (droite)
- M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
- scale.*[ sin(α), sin(α)]'], O);
- satline!(pl, M[1,:], M[2,:], subplot=subplot) # haut
- M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
- -scale.*[ sin(α), sin(α)]'], O);
- satline!(pl, M[1,:], M[2,:], subplot=subplot) # bas
-
- # Panneau
- ep = (lb-3*lp)/6
- panneau = [0 Lp Lp 0
- 0 0 lp lp]
-
- ey = 3*eb # eloignement des panneaux au bras
- vy = [cos(β-π/2); sin(β-π/2)] .* ey
- v0 = [0; yb]
- v1 = 2*ep*[cos(β); sin(β)]
- v2 = (3*ep+lp)*[cos(β); sin(β)]
- v3 = (4*ep+2*lp)*[cos(β); sin(β)]
-
- pa1 = T(R(β-π/2)*panneau, v0+v1+vy); pa = T(Rθ*pa1, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa2 = T(R(β-π/2)*panneau, v0+v2+vy); pa = T(Rθ*pa2, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa3 = T(R(β-π/2)*panneau, v0+v3+vy); pa = T(Rθ*pa3, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
-
- pa4 = SA(pa1, β, [xb; yb]); pa = T(Rθ*pa4, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa5 = SA(pa2, β, [xb; yb]); pa = T(Rθ*pa5, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa6 = SA(pa3, β, [xb; yb]); pa = T(Rθ*pa6, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
-
- pa7 = SA(pa1, 0, [0; 0]); pa = T(Rθ*pa7, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa8 = SA(pa2, 0, [0; 0]); pa = T(Rθ*pa8, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa9 = SA(pa3, 0, [0; 0]); pa = T(Rθ*pa9, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
-
- pa10 = SA(pa7, -β, [xb; -yb]); pa = T(Rθ*pa10, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa11 = SA(pa8, -β, [xb; -yb]); pa = T(Rθ*pa11, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
- pa12 = SA(pa9, -β, [xb; -yb]); pa = T(Rθ*pa12, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
-
- # Dessin flamme
- if thrust
- lthrust_max = 8000.0
- origin = T(Rθ*[Δx, 0.0], O)
- thrust_ = thrust_val*Rθ*[lthrust_max, 0.0]
- quiver!(pl, [origin[1]], [origin[2]], color=:red,
- quiver=([thrust_[1]], [thrust_[2]]), linewidth=1, subplot=subplot)
- end
-
- end;
-
- # --------------------------------------------------------------------------------------------
- #
- # Stars
- #
- # --------------------------------------------------------------------------------------------
- @userplot StarsPlot
- @recipe function f(cp::StarsPlot)
- xo, yo, Δx, Δy, I, A, m, n = cp.args
- seriestype --> :scatter
- seriescolor --> :white
- seriesalpha --> A
- markerstrokewidth --> 0
- markersize --> 2*I[3]
- xo.+I[2].*Δx./n, yo.+I[1].*Δy./m
- end
-
- mutable struct Stars
- s
- i
- m
- n
- a
- end
-
- # Etoiles
- function stars(ρ)
- n = 200 # nombre de colonnes
- m = int(ρ*n) # nombre de lignes
- d = 0.03
- s = sprand(m, n, d)
- i = findnz(s)
- s = dropzeros(s)
- # alpha
- amin = 0.6
- amax = 1.0
- Δa = amax-amin
- a = amin.+rand(length(i[3])).*Δa
- return Stars(s, i, m, n, a)
- end;
-
- # --------------------------------------------------------------------------------------------
- #
- # Trajectory
- #
- # --------------------------------------------------------------------------------------------
- @userplot TrajectoryPlot
- @recipe function f(cp::TrajectoryPlot)
- t, x, y, tf = cp.args
- n = argmin(abs.(t.-tf))
- inds = 1:n
- seriescolor --> :white
- linewidth --> range(0.0, 3.0, length=n)
- seriesalpha --> range(0.0, 1.0, length=n)
- aspect_ratio --> 1
- label --> false
- x[inds], y[inds]
- end
-
- function trajectoire!(px, times, x1, x2, θ, F, thrust, t_current, cx, cy, pt)
-
- i_current = argmin(abs.(times.-t_current))
-
- # Trajectoire
- if t_current>0 && pt[i_current]
- trajectoryplot!(px, times, cx.+x1, cy.+x2, t_current)
- end
-
- # Satellite
- satellite!(px, position=[cx+x1[i_current]; cy+x2[i_current]], scale=sc_sat,
- rotate=θ[i_current], thrust=thrust[i_current], thrust_val=F[i_current])
-
- end;
-
- # --------------------------------------------------------------------------------------------
- #
- # Background
- #
- # --------------------------------------------------------------------------------------------
- # Décor
- function background(w, h, xmin, xmax, ymin, ymax,
- X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
- cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, opi, opf, times, time_current)
-
- # Fond
- px = plot(background_color=:gray8, legend = false, aspect_ratio=:equal,
- size = (w, h), framestyle = :none,
- left_margin=-25mm, bottom_margin=-10mm, top_margin=-10mm, right_margin=-10mm,
- xlims=(xmin, xmax), ylims=(ymin, ymax) #,
- #camera=(0, 90), xlabel = "x", ylabel= "y", zlabel = "z", right_margin=25mm
- )
-
- # Etoiles
- starsplot!(px, xmin, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
-
- starsplot!(px, xmin-Δx, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin-Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin+Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin-Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
- starsplot!(px, xmin+Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
-
- # Orbite initiale
- nn = length(X2_orb_init)
- plot!(px, cx.+X1_orb_init, cy.+X2_orb_init, color = :olivedrab1, linewidth=2, alpha=opi)
-
- # Orbite finale
- nn = length(X2_orb_init)
- plot!(px, cx.+X1_orb_arr, cy.+X2_orb_arr, color = :turquoise1, linewidth=2, alpha=opf)
-
- # Terre
- θ = range(0.0, 2π, length=100)
- rT = rayon_Terre #- 1000
- xT = cx .+ rT .* cos.(θ)
- yT = cy .+ rT .* sin.(θ)
- nn = length(xT)
- plot!(px, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
-
- # Soleil
- i_current = argmin(abs.(times.-time_current))
- e = π/6
- β = range(π/4+e, π/4-e, length=length(times))
- ρ = sqrt((0.8*xmax-cx)^2+(0.8*ymax-cy)^2)
- cxS = cx + ρ * cos(β[i_current])
- cyS = cy + ρ * sin(β[i_current])
- θ = range(0.0, 2π, length=100)
- rS = 2000
- xS = cxS .+ rS .* cos.(θ)
- yS = cyS .+ rS .* sin.(θ)
- plot!(px, xS, yS, color = :gold, seriestype=:shape, linewidth=0)
-
- # Point de départ
- #plot!(px, [cx+x0[1]], [cy+x0[2]], seriestype=:scatter, color = :white, markerstrokewidth=0)
-
- # Cadre
- #plot!(px, [xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], color=:white)
-
- # Angle
- #plot!(px, camera=(30,90))
-
- return px
-
- end;
-
- function panneau_control!(px, time_current, times, u, F_max, subplot_current)
-
- # We set the (optional) position relative to top-left of the 1st subplot.
- # The call is `bbox(x, y, width, height, origin...)`,
- # where numbers are treated as "percent of parent".
- Δt_control = 2 #0.1*times[end]
-
- xx = 0.08
- yy = 0.1
- plot!(px, times, F_max.*u[1,:],
- inset = (1, bbox(xx, yy, 0.2, 0.15, :top, :left)),
- subplot = subplot_current,
- bg_inside = nothing,
- color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
- xguidefontsize=14, legendfontsize=14,
- ylims=(-F_max,F_max),
- xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
- ylabel="u₁ [N]", label=""
- )
-
- plot!(px,
- [time_current, time_current], [-F_max,F_max],
- subplot = subplot_current, label="")
-
- plot!(px, times, F_max.*u[2,:],
- inset = (1, bbox(xx, yy+0.2, 0.2, 0.15, :top, :left)),
- #ticks = nothing,
- subplot = subplot_current+1,
- bg_inside = nothing,
- color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
- xguidefontsize=14, legendfontsize=14,
- ylims=(-F_max,F_max),
- xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
- ylabel="u₂ [N]", xlabel="temps [h]", label=""
- )
-
- plot!(px,
- [time_current, time_current], [-F_max,F_max],
- subplot = subplot_current+1, label="")
-
- return subplot_current+2
- end;
-
- @enum PB tfmin=1 consomin=2
-
- function panneau_information!(px, subplot_current, time_current, i_current,
- x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr;
- pb=tfmin, tf_min=0.0, conso=[])
-
- # panneaux information
- xx = 0.06
- yy = 0.3
- plot!(px,
- inset = (1, bbox(xx, yy, 0.2, 0.15, :bottom, :left)),
- subplot=subplot_current,
- bg_inside = nothing, legend=:false,
- framestyle=:none
- )
-
- a = 0.0
- b = 0.0
- Δtxt = 0.3
-
- txt = @sprintf("Temps depuis départ [h] = %3.2f", time_current)
- plot!(px, subplot = subplot_current,
- annotation=((a, b+3*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- txt = @sprintf("Vitesse satellite [km/h] = %5.2f", sqrt(v1[i_current]^2+v2[i_current]^2))
- plot!(px, subplot = subplot_current,
- annotation=((a, b+2*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- txt = @sprintf("Distance à la Terre [km] = %5.2f", sqrt(x1[i_current]^2+x2[i_current]^2)-rayon_Terre)
- plot!(px, subplot = subplot_current,
- annotation=((a, b+Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- txt = @sprintf("Poussée maximale [N] = %5.2f", F_max)
- plot!(px, subplot = subplot_current,
- annotation=((a, b-0*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- txt = @sprintf("Durée du transfert [h] = %5.2f", tf_transfert)
- plot!(px, subplot = subplot_current,
- annotation=((a, b-1*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- if pb==consomin
- txt = @sprintf("Durée et consommation comparées")
- plot!(px, subplot = subplot_current,
- annotation=((a, b-2*Δtxt, txt)),
- annotationcolor=:green, annotationfontsize=16,
- annotationhalign=:left)
-
- txt = @sprintf("au problème en temps minimal :")
- plot!(px, subplot = subplot_current,
- annotation=((a, b-2.75*Δtxt, txt)),
- annotationcolor=:green, annotationfontsize=16,
- annotationhalign=:left)
-
- txt = @sprintf("Durée du transfert [relatif] = %5.2f", tf_transfert/tf_min)
- plot!(px, subplot = subplot_current,
- annotation=((a, b-3.75*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
-
- txt = @sprintf("Consommation [relative] = %5.2f", conso[i_current]./tf_min)
- plot!(px, subplot = subplot_current,
- annotation=((a, b-4.75*Δtxt, txt)),
- annotationcolor=:white, annotationfontsize=14,
- annotationhalign=:left)
- end
-
- subplot_current = subplot_current+1
- plot!(px,
- inset = (1, bbox(0.3, 0.03, 0.5, 0.1, :top, :left)),
- subplot = subplot_current,
- bg_inside = nothing, legend=:false, aspect_ratio=:equal,
- xlims=(-4, 4),
- framestyle=:none
- )
-
- rmax = 3*maximum(sqrt.(X1_orb_init.^2+X2_orb_init.^2))
- plot!(px, subplot = subplot_current,
- 0.04.+X1_orb_init./rmax, X2_orb_init./rmax.-0.03,
- color = :olivedrab1, linewidth=2)
-
- rmax = 7*maximum(sqrt.(X1_orb_arr.^2+X2_orb_arr.^2))
- plot!(px, subplot = subplot_current,
- 3.3.+X1_orb_arr./rmax, X2_orb_arr./rmax.-0.03,
- color = :turquoise1, linewidth=2)
-
- satellite!(px, subplot=subplot_current,
- position=[0.67; 0.28], scale=0.08,
- rotate=π/2)
-
- # Terre
- θ = range(0.0, 2π, length=100)
- rT = 0.1
- xT = 3.55 .+ rT .* cos.(θ)
- yT = 0.28 .+ rT .* sin.(θ)
- plot!(px, subplot=subplot_current, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
-
- s1 = "Transfert orbital du satellite autour de la Terre\n"
- s2 = "de l'orbite elliptique vers l'orbite circulaire\n"
- if pb==tfmin
- s3 = "en temps minimal"
- elseif pb==consomin
- s3 = "en consommation minimale"
- else
- error("pb pb")
- end
- txt = s1 * s2 * s3
- plot!(px, subplot = subplot_current,
- annotation=((0, 0, txt)),
- annotationcolor=:white, annotationfontsize=18,
- annotationhalign=:center)
-
- # Realisation
- subplot_current = subplot_current+1
- xx = 0.0
- yy = 0.02
- plot!(px,
- inset = (1, bbox(xx, yy, 0.12, 0.05, :bottom, :right)),
- subplot = subplot_current,
- bg_inside = nothing, legend=:false,
- framestyle=:none
- )
-
- s1 = "Réalisation : Olivier Cots (2022)"
- txt = s1
- plot!(px, subplot = subplot_current,
- annotation=((0, 0, txt)),
- annotationcolor=:gray, annotationfontsize=8, annotationhalign=:left)
-
- return subplot_current+1
-
- end;
-
-end # module
\ No newline at end of file