foo

c18abfc0 · Olivier Cots · c8f25ee7 · c18abfc0
Commit c18abfc0 authored 2 years ago by Olivier Cots
--- a/tp/direct-indirect-part1.ipynb
+++ b/tp/direct-indirect-part1.ipynb
@@ -41,7 +41,7 @@
    "    \\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",
+    "        x(0) = 1, \\quad x(t_f) = 1/2\n",
    "    \\end{array}\n",
    "    \\right. \n",
    "$$\n",

 %% Cell type:markdown id: tags:

 [<img src="https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/inp-enseeiht.jpg" alt="N7" height="100"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)

 # Méthode directe et tir simple structuré

 - Date : 2023
 - Durée approximative : 1.5 séance

 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 id: tags:

 ``` julia
 using JuMP, Ipopt # pour la méthode directe
 using DifferentialEquations, NLsolve, ForwardDiff # pour la méthode indirecte
 using Plots # pour les graphiques
 include("utils.jl"); # fonctions utilitaires
 ```

 %% Cell type:markdown id: tags:

 ## Introduction

 On considère le problème de contrôle optimal suivant :

 $$
    \left\{
    \begin{array}{l}
        \displaystyle \min_{x, u} \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \\[1.0em]
        \dot{x}(t) = \displaystyle u(t), \quad  |u(t)| \le 1, \quad t \in [0, t_f] \text{ p.p.},    \\[1.0em]
-        x(0) = 1, \quad x(t_f) = 1/2.
+        x(0) = 1, \quad x(t_f) = 1/2
    \end{array}
    \right.
 $$
 où $t_f$ = 2.

 ✏️ **Exercice 1.**

 1. Appliquez le PMP. Que pouvez-vous dire du contrôle maximisant ?
 2. Peut-on appliquer simplement la méthode de tir simple vu au TP précédent ?

 %% Cell type:markdown id: tags:

 ## Méthode directe

 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.

 ✏️ **Exercice 2.**

 - 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)$.

 %% Cell type:markdown id: tags:

 **Description de la méthode directe.**

 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.

 Pour cela :

 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.
 2. On discrétise l'état, le contrôle et le coût sur cette grille et on note
 $$
 X = \{ (x_i, u_i, c_i) ~|~ i \in \{1, \dots, N+1\}\}
 $$
 l'ensemble des variables d'optimisation du problème discrétisé.

 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
 $$
 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\}.
 $$

 ✏️ ️**Exercice 3.** Définissez pour le problème discrétisé :

 1. l'objectif à minimiser en fonction d'une ou plusieurs composantes de $X$ bien choisies.
 2. les contraintes d'inégalités associées au contrôle.
 3. les contraintes initiales et finales associées à la variable d'état.
 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)).

 %% Cell type:markdown id: tags:

 ✏️ ️**Exercice 4.**

 - Modifiez le code suivant afin de résoudre le problème par la méthode directe que vous venez de décrire.

 **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 id: tags:

 ``` julia
 # Création du modèle JuMP, Utilisation de Ipopt comme solver
 sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5))
 set_optimizer_attribute(sys,"tol",1e-8)
 set_optimizer_attribute(sys,"constr_viol_tol",1e-6)
 set_optimizer_attribute(sys,"max_iter",1000)

 #######
 ####### DEBUT A MODIFIER
 #######
 ####### Les ... sont a remplacer !
 #######
 ####### Attention, on ecrit le probleme sous la forme d'un probleme de Mayer
 #######

 # Paramètres
 t0 = 0    # temps initial
 tf = 0    # temps final
 c0 = 0    # coût initial
 x0 = 0    # état initial
 xf = 0   # état final

 N  = 50    # taille de la grille
 Δt = (tf-t0)/N  # pas de temps

 @variables(sys, begin
    ...                # coût
    ...                # état
    ...       # control
 end)

 # Objectif
 @objective(sys, Min, ...)

 # Conditions initiales et finales
 @constraints(sys, begin
    con_c0, ... # contraint sur le coût initial
    con_x0, ... # contraint sur l'état initial
    con_xf, ... # contraint sur l'état final
 end)

 # Contraintes de dynamique avec le schéma de Crank-Nicolson
 @NLconstraints(sys, begin
    con_dc[j=1:N], ... # contraint différentielle sur le coût
    con_dx[j=1:N], ... # contraint différentielle sur l'état
 end);


 #######
 ####### FIN A MODIFIER
 #######
 ```

 %% Cell type:code id: tags:

 ``` julia
 # Solve for the control and state
 println("Solving...")
 optimize!(sys)
 println()

 # Display results
 if termination_status(sys) == MOI.OPTIMAL
    println("  Solution is optimal")
 elseif  termination_status(sys) == MOI.LOCALLY_SOLVED
    println("  (Local) solution found")
 elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys)
    println("  Solution is suboptimal due to a time limit, but a primal solution is available")
 else
    error("  The model was not solved correctly.")
 end
 println("  objective value = ", objective_value(sys))
 println()

 # Retrieves values (including duals)
 c  = value.(c)[:]
 x  = value.(x)[:]
 u  = value.(u)[:]
 t  = (0:N) * value.(Δt)

 pc0 = dual(con_c0)
 px0 = dual(con_x0)
 pxf = dual(con_xf)

 if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end
 if(px0*dual(con_dx[1])<0); px0 = -px0; end
 if(pxf*dual(con_dx[N])<0); pxf = -pxf; end

 if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle
    sign = -1.0
 else
    sign =  1.0
 end

 pc = [ dual(con_dc[i]) for i in 1:N ]
 px = [ dual(con_dx[i]) for i in 1:N ]

 pc = sign * [pc0; pc[1:N]]
 px = sign * [px0; (px[1:N-1]+px[2:N])/2; pxf]; # We add the multiplier from the limit conditions
 ```

 %% Cell type:code id: tags:

 ``` julia
 x_plot = plot(t, x, xlabel = "t", ylabel = "x", legend = false)
 u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,400), linetype=:steppre)
 px_plot = plot(t, px, xlabel = "t", ylabel = "p", legend = false)
 display(plot(x_plot, px_plot, layout = (1,2), size=(800,400)))
 display(u_plot)
 ```

 %% Cell type:markdown id: tags:

 ✏️ ️**Exercice 5.**

 1. Commentez les résultats.
 2. Modifiez la tolérance de l'optimiseur ainsi que le nombre de points de discrétisation. Commentaires.
 3. Décrivez la structure du contrôle optimal en fonction du temps.