From a6087c4c36afec6da85726e2edacd187f494a3bc Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 8 May 2024 15:52:43 +0200
Subject: [PATCH] add TP4
---
README.md | 4 +-
tp/direct-indirect.ipynb | 430 +++++++++++++++++++++++++++++++++++++++
2 files changed, 432 insertions(+), 2 deletions(-)
create mode 100644 tp/direct-indirect.ipynb
diff --git a/README.md b/README.md
index a23d3d8..754b46e 100644
--- a/README.md
+++ b/README.md
@@ -12,8 +12,8 @@ Cours de contrôle optimal de l'ENSEEIHT pour le département Sciences du Numér
* [Sujet 1 : tp/ode.ipynb](tp/ode.ipynb) - Intégration numérique et systèmes dynamiques
* [Sujet 2 : tp/simple-shooting.ipynb](tp/simple-shooting.ipynb) - Tir simple indirect
-* Sujet 3 : en cours de réalisation
-* [Sujet 4 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul numérique de dérivées
+* [Sujet 3 : tp/derivative.ipynb](tp/derivative.ipynb) - Calcul numérique de dérivées
+* [Sujet 4 : tp/direct-indirect.ipynb](tp/direct-indirect.ipynb) - Méthode directe et tir avec structure sur le contrôle
* Sujet 5 : plus tard
**Notes de cours supplémentaires - ressources**
diff --git a/tp/direct-indirect.ipynb b/tp/direct-indirect.ipynb
new file mode 100644
index 0000000..3e6a38e
--- /dev/null
+++ b/tp/direct-indirect.ipynb
@@ -0,0 +1,430 @@
+{
+ "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
+}
--
GitLab