Compare revisions

Olivier Cots · Olivier Cots · 73562ae3 · 73562ae3 · fb80aa24 · fb80aa24
--- a/README.md
+++ b/README.md
@@ -6,7 +6,7 @@ Cours de contrôle optimal de l'ENSEEIHT pour le département Sciences du Numér
 **CTD**
-* [Notes de cours - contrôle optimal](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/notes-co-2024.pdf)
+* [Notes de cours - contrôle optimal](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/notes-co-2025.pdf)
 **TP**
@@ -24,10 +24,6 @@ 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
-Projet
-* [Projet : tp/transfert-orbital.ipynb](tp/transfert-orbital.ipynb) - Transfert orbital
 **Notes de cours supplémentaires - ressources**
 * [Automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf)

--- a/notes-co-2024.pdf
+++ b/notes-co-2024.pdf
--- a/stages/.gitkeep
+++ b/stages/.gitkeep
--- a/stages/2024-stage-M1-cots.pdf
+++ b/stages/2024-stage-M1-cots.pdf
--- a/stages/Internship_subject_2.pdf
+++ b/stages/Internship_subject_2.pdf
--- a/tp/transfert-orbital.ipynb
+++ b/tp/transfert-orbital.ipynb
-%% Cell type:markdown id: tags:
-[<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)
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#FAA299;
-            border:2px solid #BF381B;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-* Nom :
-* Prénom :
-</div>
-%% Cell type:markdown id: tags:
-# Transfert orbital
-Ce sujet (ou TP-projet) est à rendre (voir la date sur moodle et les modalités du rendu) et sera évalué pour faire partie de la note finale de la matière Contrôle Optimal.
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#afa;
-            border:2px solid #bbffbb;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1.5em;
-            text-align:center;">
-Transfert orbital à temps minimal
-</div>
-On considère le problème du transfert d'un satellite d'une orbite initiale à l'orbite géostationnaire à temps minimal. Ce problème s'écrit comme un problème de contrôle optimal sous la forme
-$$
-\left\lbrace
-\begin{array}{l}
-    \min J(x, u, t_f) = t_f \\[1.0em]
-    \ \ \dot{x}_{1}(t) = ~ x_{3}(t)  \\[0.5em]
-    \ \ \dot{x}_{2}(t) = ~ x_{4}(t)  \\[0.5em]
-    \ \ \dot{x}_{3}(t) =  -\dfrac{\mu\, x_{1}(t)}{r^{3}(x(t))} + u_{1}(t)  \\[1em]
-    \ \ \dot{x}_{4}(t) =  -\dfrac{\mu\, x_{2}(t)}{r^{3}(x(t))} + u_{2}(t), ~~ ||u(t)|| \leq \gamma_\mathrm{max}, ~~ t \in [0,t_f] ~~ \text{p.p.}, ~~ u(t) = (u_1(t),u_2(t)),  \\[1.5em]
-    \ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\[1em]
-    \ \ r(x(t_f)) = r_{f}, ~~ x_{3}(t_f)=-\sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4}           (t_f)= \sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{1}(t_f), \\
-\end{array}
-\right.
-$$
-avec $r(x)=\sqrt{x_{1}^{2}+x_{2}^{2}}$.
-Les unités choisies sont le kilomètre pour les distances et l'heure pour les  temps. On donne les paramètres suivants :
-$$
-\mu=5.1658620912 \times 10^{12} \ \mathrm{km}^{3}.\mathrm{h}^{-2}, \quad r_{f} = 42165 \ \mathrm{km}.
-$$
-Le paramètre $\gamma_\mathrm{max}$ dépend de la poussée maximale $F_\mathrm{max}$ suivant la relation :
-$$
-\gamma_\mathrm{max} = \frac{F_\mathrm{max}\times 3600^2}{m}
-$$
-où m est la masse du satellite qu'on fixe à $m=2000\ \mathrm{kg}$.
-%% Cell type:markdown id: tags:
-### Résolution via du tir simple indirect
-%% Cell type:code id: tags:
-``` julia
-using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra
-include("utils.jl"); # fonctions utilitaires
-```
-%% Cell type:code id: tags:
-``` julia
-# Les constantes du pb
-x0     = [-42272.67, 0, 0, -5796.72] # état initial
-μ      = 5.1658620912*1e12
-rf     = 42165
-F_max  = 100
-γ_max  = F_max*3600^2/(2000*10^3)
-t0     = 0
-rf3    = rf^3
-α      = sqrt(μ/rf3);
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 1:_**
-1. Donner le pseudo-hamiltonien associé au problème de contrôle optimal.
-2. Donner le pseudo système hamiltonien
-$$
-    \vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-    -\frac{\partial H}{\partial x}(x, p, u) \right).
-$$
-3. Calculer le contrôle maximisant. On supposera que $(p_3, p_4)\neq (0,0)$.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-✏️ **_Question 2:_** Compléter le code suivant contenant le pseudo-hamiltonien, le pseudo système hamiltonien et le contrôle maximisant.
-%% Cell type:code id: tags:
-``` julia
-#####
-##### A COMPLETER
-# pseudo-Hamiltonien
-function H(x, p, u)
-    h = 0
-    return h
-end
-# pseudo système hamiltonien
-function Hv(x, p, u)
-    n     = size(x, 1)
-    dx    = zeros(eltype(x), n)
-    dp    = zeros(eltype(x), n)
-    return dx, dp
-end
-# Contrôle maximisant
-function control(p)
-    u    = zeros(eltype(p),2)
-    return u
-end
-#####
-##### FIN A COMPLETER
-# flot hamiltonien pour le calcul des extrémales
-f = Flow((x, p) -> Hv(x, p, control(p)));
-```
-%% Cell type:markdown id: tags:
-On note
-$
-    \alpha := \sqrt{\frac{\mu}{r_f^3}}.
-$
-La condition terminale peut se mettre sous la forme $c(x(t_f)) = 0$, avec $c \colon \mathbb{R}^4 \to \mathbb{R}^3$.
-✏️ **_Question 3:_** Donner l'expression de $c(x)$.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-Le temps final étant libre, on a la condition au temps final
-$$
-    H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(on se place dans le cas normal)}
-$$
-De plus, la condition de transversalité
-$$
-p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
-$$
-conduit à la relation suivante (où $\lambda$ n'apparaît plus)
-$$
-\Phi(x(t_f), p(t_f)) := x_2(t_f) \Big( p_1(t_f) + \alpha\, p_4(t_f) \Big) - x_1(t_f) \Big( p_2(t_f) - \alpha\, p_3(t_f) \Big) = 0.
-$$
-En considérant la condition aux limites, la condition finale sur le pseudo-hamiltonien et la condition de transversalité, la fonction de tir simple est donnée par
-\begin{equation*}
-    \begin{array}{rlll}
-        S \colon    & \mathbb{R}^5          & \longrightarrow   & \mathbb{R}^5 \\
-        & (p_0, t_f)      & \longmapsto       &
-        S(p_0, t_f) := \begin{pmatrix}
-            c(x(t_f, x_0, p_0)) \\[0.5em]
-            \Phi(z(t_f, x_0, p_0)) \\[0.5em]
-            H(z(t_f, x_0, p_0), u(z(t_f, x_0, p_0))) - 1
-        \end{pmatrix}
-    \end{array}
-\end{equation*}
-où $z(t_f, x_0, p_0)$ est la solution au temps $t_f$ du pseudo système hamiltonien bouclé par
-le contrôle maximisant, partant au temps $t_0=0$ du point $(x_0, p_0)$. On rappelle que l'on note
-$z=(x, p)$.
-%% Cell type:markdown id: tags:
-✏️ **_Question 4:_** Compléter le code suivant de la fonction de tir.
-%% Cell type:code id: tags:
-``` julia
-#####
-##### A COMPLETER
-# Fonction de tir
-function shoot(p0, tf)
-    s = zeros(eltype(p0), 5)
-    ...
-    return s
-end;
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 5:_** Expliquer simplement ce que fait le code suivant.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:code id: tags:
-``` julia
-# Itéré initial
-y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4]   # pour F_max = 100N
-# Jacobienne de la fonction de tir
-foo(y)  = shoot(y[1:4], y[5])
-jfoo(y) = ForwardDiff.jacobian(foo, y)
-# Résolution de shoot(p0, tf) = 0
-nl_sol = nlsolve(foo, jfoo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true);
-# On récupère la solution si convergence
-if converged(nl_sol)
-    p0_sol_100 = nl_sol.zero[1:4]
-    tf_sol_100 = nl_sol.zero[5]
-    println("\nFinal solution:\n", nl_sol.zero)
-else
-    error("Not converged")
-end
-```
-%% Cell type:code id: tags:
-``` julia
-# Fonction d'affichage d'une solution
-function plot_solution(p0, tf)
-    # On trace l'orbite de départ et d'arrivée
-    gr(dpi=300, size=(500,400), thickness_scaling=1)
-    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   = 101
-    Theta     = range(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)
-    # Orbite initiale
-    for  i in 1:n_theta
-        theta = Theta[i]
-        r_orb = p_orb/(1+e_ellipse*cos(theta));
-        X1_orb_init[i] = r_orb*cos(theta);
-        X2_orb_init[i] = r_orb*sin(theta);
-    end
-    # Orbite d'arrivée
-    for  i in 1:n_theta
-        theta = Theta[i]
-        X1_orb_arr[i] = rf*cos(theta) ;
-        X2_orb_arr[i] = rf*sin(theta);
-    end;
-    # Calcul de la trajectoire
-    ode_sol  = f((t0, tf), x0, p0)
-    t  = ode_sol.t
-    x1 = [ode_sol[1, j] for j in 1:size(t, 1) ]
-    x2 = [ode_sol[2, j] for j in 1:size(t, 1) ]
-    u  = zeros(2, length(t))
-    for j in 1:size(t, 1)
-        u[:,j] = control(ode_sol[5:8, j])
-    end
-    px = plot(x1, x2, color=:red, legend=:best, label="Trajectoire",
-        border=:none, size = (800,400), aspect_ratio=:equal)
-    plot!(px, X1_orb_init, X2_orb_init, color=:green, label="Orbite initiale")
-    plot!(px, X1_orb_arr, X2_orb_arr, color=:blue, label="Orbite d'arrivée")
-    plot!(px, [x0[1]], [x0[2]], seriestype=:scatter, color=:green, markersize = 5, label="Départ")
-    xf = ode_sol[1:2, end]
-    plot!(px, [xf[1]], [xf[2]], seriestype=:scatter, color=:red, markersize = 5, label="Arrivée")
-    plot!(px, [0.0], [0.0], color = :blue, seriestype=:scatter, markersize = 25, label="Terre")
-    pu1 = plot(t, u[1,:], color=:red, label="u₁(t)", legend=:best)
-    pu2 = plot(t, u[2,:], color=:red, label="u₂(t)", legend=:best)
-    display(plot(pu1, pu2, layout = (1,2), size = (800,400)))
-    display(px)
-end;
-```
-%% Cell type:code id: tags:
-``` julia
-# On affiche la solution pour Fmax = 100
-plot_solution(p0_sol_100, tf_sol_100)
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 6:_** Que vaut le temps final $t_f$ (à 5 digits près) ? Combien de révolutions complètes autour de la Terre a réalisé le satellite.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-### Etude du temps de transfert en fonction de la poussée maximale
-✏️ **_Question 7:_**
-* A l'aide de ce que vous avez fait précédemment, déterminer $t_f$ (attention, penser à bien itinialiser la
-méthode de tir) pour
-$$
-    F_\mathrm{max} \in \{100, 90, 80, 70, 60, 50, 40, 30, 20 \}.
-$$
-* Tracer ensuite $t_f$ en fonction de $F_\mathrm{max}$ et commenter le résultat.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte et donner votre commentaire)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:code id: tags:
-``` julia
-# Les différentes valeurs de poussées
-F_max_span = range(100, stop=20, length=11);
-γ_max_span = [F_max_span[j]*3600^2/(2000*10^3) for j in 1:size(F_max_span,1)];
-```
-%% Cell type:code id: tags:
-``` julia
-# Solution calculée précédemment
-y_guess = [p0_sol_100; [tf_sol_100]]
-# Pour le stockage des solutions
-tf_sols = zeros(length(γ_max_span))     # vecteur des temps de transfert
-p0_sols = zeros(4, length(γ_max_span))  # matrice des co-états initiaux
-#####
-##### A COMPLETER
-...
-```
-%% Cell type:code id: tags:
-``` julia
-# Affichage de tf en fonction de la poussée maximale
-plot(F_max_span, tf_sols, aspect_ratio=:equal, legend=:best, label="tf(Fmax)")
-```
-%% Cell type:code id: tags:
-``` julia
-# Plots sol pour F_max = 20N
-i = 11
-γ_max = γ_max_span[i]
-plot_solution(p0_sols[:, i], tf_sols[i])
-```
-%% Cell type:markdown id: tags:
-### Animation
-Juste pour s'amuser
-%% Cell type:code id: tags:
-``` julia
-include("space.jl"); using .Space
-```
-%% Cell type:code id: tags:
-``` julia
-# solution choice
-i = 11
-γ_max = γ_max_span[i] # γ_max must be updated for the use of the flow
-# animation
-# nFrame = 100; d = 30; fps = floor(Int, nFrame/d) => 2 minutes and 30 seconds of computation
-nFrame = 100; d = 30; fps = floor(Int, nFrame/d)
-Space.animation(p0_sols[:, i], tf_sols[i], f, γ_max, nFrame=nFrame, fps=fps)
-```
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#afa;
-            border:2px solid #bbffbb;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1.5em;
-            text-align:center;">
-Transfert orbital à consommation minimale
-</div>
-%% Cell type:markdown id: tags:
-On considère maintenant le problème du transfert d'un satellite d'une orbite initiale à une orbite géostationnaire à temps final fixé en cherchant à minimiser la consommation de carburant. Ce problème s'écrit comme un problème de contrôle optimal sous la forme
-$$
-\left\lbrace
-\begin{array}{l}
-    \displaystyle \min J(x, u) = \int_{0}^{t_f} \Vert u(t)\Vert \mathrm{d}t \\[1.0em]
-    \ \ \dot{x}_{1}(t) = ~ x_{3}(t)  \\[0.5em]
-    \ \ \dot{x}_{2}(t) = ~ x_{4}(t)  \\[0.5em]
-    \ \ \dot{x}_{3}(t) =  -\dfrac{\mu\, x_{1}(t)}{r^{3}(x(t))} + u_{1}(t)  \\[1em]
-    \ \ \dot{x}_{4}(t) =  -\dfrac{\mu\, x_{2}(t)}{r^{3}(x(t))} + u_{2}(t), ~~ ||u(t)|| \leq \gamma_\mathrm{max}, ~~ t \in [0,t_f] ~~ \text{p.p.}, ~~ u(t) = (u_1(t),u_2(t)),  \\[1.5em]
-    \ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\[1em]
-    \ \ r(x(t_f)) = r_{f}, ~~ x_{3}(t_f)=-\sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4}           (t_f)= \sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{1}(t_f). \\
-\end{array}
-\right.
-$$
-Toutes les constantes sont identiques au problème précédent. On considéra le problème avec $F_\mathrm{max} = 100$ et on prendre comme temps final : $t_f = 1.5 T^{100}_{\min}$, où $T^{100}_{\min}$ correspond au temps minimal solution du problème précédent pour $F_\mathrm{max} = 100$.
-**Remarques importantes.**
-* On ne considèrera que des extrémales normales, c'est-à-dire $p^0 = -1$.
-* Le problème a de la structure, le contrôle optimal est composé d'arcs où la poussée est nulle $u(t) = 0$ et où la poussée est maximale $\Vert u(t)\Vert = \gamma_\mathrm{max}$.
-* Vous devrez dans un premier temps résoudre le problème avec le coût suivant ce qui permet de le régulariser et ainsi d'utiliser une méthode de tir simple :
-$$
-\int_{0}^{t_f} \left( \Vert u(t)\Vert - \varepsilon \left( \ln(\Vert u(t)\Vert) + \ln(\gamma_\mathrm{max}-\Vert u(t)\Vert) \right) \right) \mathrm{d}t.
-$$
-%% Cell type:markdown id: tags:
-✏️ **_Question 8:_** Résoudre via du tir simple le problème régularisé pour différentes valeurs de $\varepsilon$ suffisamment petites pour déterminer la structure optimale, en commençant avec $\varepsilon$ suffisamment grand pour faire converger la méthode de tir plus facilement.
-%% Cell type:code id: tags:
-``` julia
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 9:_** Résoudre via du tir multiple le problème de transfert orbital à consommation minimale. Vous utilisez la structure obtenue précédemment et vous donnerez au solveur une bonne condition initiale à l'aide de la question précédente.
-%% Cell type:code id: tags:
-``` julia
-```
-%% Cell type:markdown id: tags:
-[<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)
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#FAA299;
-            border:2px solid #BF381B;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-* Nom :
-* Prénom :
-</div>
-%% Cell type:markdown id: tags:
-# Transfert orbital
-Ce sujet (ou TP-projet) est à rendre (voir la date sur moodle et les modalités du rendu) et sera évalué pour faire partie de la note finale de la matière Contrôle Optimal.
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#afa;
-            border:2px solid #bbffbb;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1.5em;
-            text-align:center;">
-Transfert orbital à temps minimal
-</div>
-On considère le problème du transfert d'un satellite d'une orbite initiale à l'orbite géostationnaire à temps minimal. Ce problème s'écrit comme un problème de contrôle optimal sous la forme
-$$
-\left\lbrace
-\begin{array}{l}
-    \min J(x, u, t_f) = t_f \\[1.0em]
-    \ \ \dot{x}_{1}(t) = ~ x_{3}(t)  \\[0.5em]
-    \ \ \dot{x}_{2}(t) = ~ x_{4}(t)  \\[0.5em]
-    \ \ \dot{x}_{3}(t) =  -\dfrac{\mu\, x_{1}(t)}{r^{3}(x(t))} + u_{1}(t)  \\[1em]
-    \ \ \dot{x}_{4}(t) =  -\dfrac{\mu\, x_{2}(t)}{r^{3}(x(t))} + u_{2}(t), ~~ ||u(t)|| \leq \gamma_\mathrm{max}, ~~ t \in [0,t_f] ~~ \text{p.p.}, ~~ u(t) = (u_1(t),u_2(t)),  \\[1.5em]
-    \ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\[1em]
-    \ \ r(x(t_f)) = r_{f}, ~~ x_{3}(t_f)=-\sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4}           (t_f)= \sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{1}(t_f), \\
-\end{array}
-\right.
-$$
-avec $r(x)=\sqrt{x_{1}^{2}+x_{2}^{2}}$.
-Les unités choisies sont le kilomètre pour les distances et l'heure pour les  temps. On donne les paramètres suivants :
-$$
-\mu=5.1658620912 \times 10^{12} \ \mathrm{km}^{3}.\mathrm{h}^{-2}, \quad r_{f} = 42165 \ \mathrm{km}.
-$$
-Le paramètre $\gamma_\mathrm{max}$ dépend de la poussée maximale $F_\mathrm{max}$ suivant la relation :
-$$
-\gamma_\mathrm{max} = \frac{F_\mathrm{max}\times 3600^2}{m}
-$$
-où m est la masse du satellite qu'on fixe à $m=2000\ \mathrm{kg}$.
-%% Cell type:markdown id: tags:
-### Résolution via du tir simple indirect
-%% Cell type:code id: tags:
-``` julia
-using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra
-include("utils.jl"); # fonctions utilitaires
-```
-%% Cell type:code id: tags:
-``` julia
-# Les constantes du pb
-x0     = [-42272.67, 0, 0, -5796.72] # état initial
-μ      = 5.1658620912*1e12
-rf     = 42165
-F_max  = 100
-γ_max  = F_max*3600^2/(2000*10^3)
-t0     = 0
-rf3    = rf^3
-α      = sqrt(μ/rf3);
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 1:_**
-1. Donner le pseudo-hamiltonien associé au problème de contrôle optimal.
-2. Donner le pseudo système hamiltonien
-$$
-    \vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-    -\frac{\partial H}{\partial x}(x, p, u) \right).
-$$
-3. Calculer le contrôle maximisant. On supposera que $(p_3, p_4)\neq (0,0)$.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-✏️ **_Question 2:_** Compléter le code suivant contenant le pseudo-hamiltonien, le pseudo système hamiltonien et le contrôle maximisant.
-%% Cell type:code id: tags:
-``` julia
-#####
-##### A COMPLETER
-# pseudo-Hamiltonien
-function H(x, p, u)
-    h = 0
-    return h
-end
-# pseudo système hamiltonien
-function Hv(x, p, u)
-    n     = size(x, 1)
-    dx    = zeros(eltype(x), n)
-    dp    = zeros(eltype(x), n)
-    return dx, dp
-end
-# Contrôle maximisant
-function control(p)
-    u    = zeros(eltype(p),2)
-    return u
-end
-#####
-##### FIN A COMPLETER
-# flot hamiltonien pour le calcul des extrémales
-f = Flow((x, p) -> Hv(x, p, control(p)));
-```
-%% Cell type:markdown id: tags:
-On note
-$
-    \alpha := \sqrt{\frac{\mu}{r_f^3}}.
-$
-La condition terminale peut se mettre sous la forme $c(x(t_f)) = 0$, avec $c \colon \mathbb{R}^4 \to \mathbb{R}^3$.
-✏️ **_Question 3:_** Donner l'expression de $c(x)$.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-Le temps final étant libre, on a la condition au temps final
-$$
-    H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(on se place dans le cas normal)}
-$$
-De plus, la condition de transversalité
-$$
-p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
-$$
-conduit à la relation suivante (où $\lambda$ n'apparaît plus)
-$$
-\Phi(x(t_f), p(t_f)) := x_2(t_f) \Big( p_1(t_f) + \alpha\, p_4(t_f) \Big) - x_1(t_f) \Big( p_2(t_f) - \alpha\, p_3(t_f) \Big) = 0.
-$$
-En considérant la condition aux limites, la condition finale sur le pseudo-hamiltonien et la condition de transversalité, la fonction de tir simple est donnée par
-\begin{equation*}
-    \begin{array}{rlll}
-        S \colon    & \mathbb{R}^5          & \longrightarrow   & \mathbb{R}^5 \\
-        & (p_0, t_f)      & \longmapsto       &
-        S(p_0, t_f) := \begin{pmatrix}
-            c(x(t_f, x_0, p_0)) \\[0.5em]
-            \Phi(z(t_f, x_0, p_0)) \\[0.5em]
-            H(z(t_f, x_0, p_0), u(z(t_f, x_0, p_0))) - 1
-        \end{pmatrix}
-    \end{array}
-\end{equation*}
-où $z(t_f, x_0, p_0)$ est la solution au temps $t_f$ du pseudo système hamiltonien bouclé par
-le contrôle maximisant, partant au temps $t_0=0$ du point $(x_0, p_0)$. On rappelle que l'on note
-$z=(x, p)$.
-%% Cell type:markdown id: tags:
-✏️ **_Question 4:_** Compléter le code suivant de la fonction de tir.
-%% Cell type:code id: tags:
-``` julia
-#####
-##### A COMPLETER
-# Fonction de tir
-function shoot(p0, tf)
-    s = zeros(eltype(p0), 5)
-    ...
-    return s
-end;
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 5:_** Expliquer simplement ce que fait le code suivant.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:code id: tags:
-``` julia
-# Itéré initial
-y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4]   # pour F_max = 100N
-# Jacobienne de la fonction de tir
-foo(y)  = shoot(y[1:4], y[5])
-jfoo(y) = ForwardDiff.jacobian(foo, y)
-# Résolution de shoot(p0, tf) = 0
-nl_sol = nlsolve(foo, jfoo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true);
-# On récupère la solution si convergence
-if converged(nl_sol)
-    p0_sol_100 = nl_sol.zero[1:4]
-    tf_sol_100 = nl_sol.zero[5]
-    println("\nFinal solution:\n", nl_sol.zero)
-else
-    error("Not converged")
-end
-```
-%% Cell type:code id: tags:
-``` julia
-# Fonction d'affichage d'une solution
-function plot_solution(p0, tf)
-    # On trace l'orbite de départ et d'arrivée
-    gr(dpi=300, size=(500,400), thickness_scaling=1)
-    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   = 101
-    Theta     = range(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)
-    # Orbite initiale
-    for  i in 1:n_theta
-        theta = Theta[i]
-        r_orb = p_orb/(1+e_ellipse*cos(theta));
-        X1_orb_init[i] = r_orb*cos(theta);
-        X2_orb_init[i] = r_orb*sin(theta);
-    end
-    # Orbite d'arrivée
-    for  i in 1:n_theta
-        theta = Theta[i]
-        X1_orb_arr[i] = rf*cos(theta) ;
-        X2_orb_arr[i] = rf*sin(theta);
-    end;
-    # Calcul de la trajectoire
-    ode_sol  = f((t0, tf), x0, p0)
-    t  = ode_sol.t
-    x1 = [ode_sol[1, j] for j in 1:size(t, 1) ]
-    x2 = [ode_sol[2, j] for j in 1:size(t, 1) ]
-    u  = zeros(2, length(t))
-    for j in 1:size(t, 1)
-        u[:,j] = control(ode_sol[5:8, j])
-    end
-    px = plot(x1, x2, color=:red, legend=:best, label="Trajectoire",
-        border=:none, size = (800,400), aspect_ratio=:equal)
-    plot!(px, X1_orb_init, X2_orb_init, color=:green, label="Orbite initiale")
-    plot!(px, X1_orb_arr, X2_orb_arr, color=:blue, label="Orbite d'arrivée")
-    plot!(px, [x0[1]], [x0[2]], seriestype=:scatter, color=:green, markersize = 5, label="Départ")
-    xf = ode_sol[1:2, end]
-    plot!(px, [xf[1]], [xf[2]], seriestype=:scatter, color=:red, markersize = 5, label="Arrivée")
-    plot!(px, [0.0], [0.0], color = :blue, seriestype=:scatter, markersize = 25, label="Terre")
-    pu1 = plot(t, u[1,:], color=:red, label="u₁(t)", legend=:best)
-    pu2 = plot(t, u[2,:], color=:red, label="u₂(t)", legend=:best)
-    display(plot(pu1, pu2, layout = (1,2), size = (800,400)))
-    display(px)
-end;
-```
-%% Cell type:code id: tags:
-``` julia
-# On affiche la solution pour Fmax = 100
-plot_solution(p0_sol_100, tf_sol_100)
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 6:_** Que vaut le temps final $t_f$ (à 5 digits près) ? Combien de révolutions complètes autour de la Terre a réalisé le satellite.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:markdown id: tags:
-### Etude du temps de transfert en fonction de la poussée maximale
-✏️ **_Question 7:_**
-* A l'aide de ce que vous avez fait précédemment, déterminer $t_f$ (attention, penser à bien itinialiser la
-méthode de tir) pour
-$$
-    F_\mathrm{max} \in \{100, 90, 80, 70, 60, 50, 40, 30, 20 \}.
-$$
-* Tracer ensuite $t_f$ en fonction de $F_\mathrm{max}$ et commenter le résultat.
-%% Cell type:markdown id: tags:
-**Réponse** (double cliquer pour modifier le texte et donner votre commentaire)
-<div style="width:90%;
-            margin:10px;
-            padding:20px;
-            background-color:#CFF3F7;
-            border:2px solid #063970;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1em;"> <!-- il faut laisser une ligne vide après celle-ci -->
-Ecrire la réponse ici
-</div>
-%% Cell type:code id: tags:
-``` julia
-# Les différentes valeurs de poussées
-F_max_span = range(100, stop=20, length=11);
-γ_max_span = [F_max_span[j]*3600^2/(2000*10^3) for j in 1:size(F_max_span,1)];
-```
-%% Cell type:code id: tags:
-``` julia
-# Solution calculée précédemment
-y_guess = [p0_sol_100; [tf_sol_100]]
-# Pour le stockage des solutions
-tf_sols = zeros(length(γ_max_span))     # vecteur des temps de transfert
-p0_sols = zeros(4, length(γ_max_span))  # matrice des co-états initiaux
-#####
-##### A COMPLETER
-...
-```
-%% Cell type:code id: tags:
-``` julia
-# Affichage de tf en fonction de la poussée maximale
-plot(F_max_span, tf_sols, aspect_ratio=:equal, legend=:best, label="tf(Fmax)")
-```
-%% Cell type:code id: tags:
-``` julia
-# Plots sol pour F_max = 20N
-i = 11
-γ_max = γ_max_span[i]
-plot_solution(p0_sols[:, i], tf_sols[i])
-```
-%% Cell type:markdown id: tags:
-### Animation
-Juste pour s'amuser
-%% Cell type:code id: tags:
-``` julia
-include("space.jl"); using .Space
-```
-%% Cell type:code id: tags:
-``` julia
-# solution choice
-i = 11
-γ_max = γ_max_span[i] # γ_max must be updated for the use of the flow
-# animation
-# nFrame = 100; d = 30; fps = floor(Int, nFrame/d) => 2 minutes and 30 seconds of computation
-nFrame = 100; d = 30; fps = floor(Int, nFrame/d)
-Space.animation(p0_sols[:, i], tf_sols[i], f, γ_max, nFrame=nFrame, fps=fps)
-```
-%% Cell type:markdown id: tags:
-<div style="width:90%;
-            margin:10px;
-            padding:8px;
-            background-color:#afa;
-            border:2px solid #bbffbb;
-            border-radius:20px;
-            font-weight:bold;
-            font-size:1.5em;
-            text-align:center;">
-Transfert orbital à consommation minimale
-</div>
-%% Cell type:markdown id: tags:
-On considère maintenant le problème du transfert d'un satellite d'une orbite initiale à une orbite géostationnaire à temps final fixé en cherchant à minimiser la consommation de carburant. Ce problème s'écrit comme un problème de contrôle optimal sous la forme
-$$
-\left\lbrace
-\begin{array}{l}
-    \displaystyle \min J(x, u) = \int_{0}^{t_f} \Vert u(t)\Vert \mathrm{d}t \\[1.0em]
-    \ \ \dot{x}_{1}(t) = ~ x_{3}(t)  \\[0.5em]
-    \ \ \dot{x}_{2}(t) = ~ x_{4}(t)  \\[0.5em]
-    \ \ \dot{x}_{3}(t) =  -\dfrac{\mu\, x_{1}(t)}{r^{3}(x(t))} + u_{1}(t)  \\[1em]
-    \ \ \dot{x}_{4}(t) =  -\dfrac{\mu\, x_{2}(t)}{r^{3}(x(t))} + u_{2}(t), ~~ ||u(t)|| \leq \gamma_\mathrm{max}, ~~ t \in [0,t_f] ~~ \text{p.p.}, ~~ u(t) = (u_1(t),u_2(t)),  \\[1.5em]
-    \ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\[1em]
-    \ \ r(x(t_f)) = r_{f}, ~~ x_{3}(t_f)=-\sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4}           (t_f)= \sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{1}(t_f). \\
-\end{array}
-\right.
-$$
-Toutes les constantes sont identiques au problème précédent. On considéra le problème avec $F_\mathrm{max} = 100$ et on prendre comme temps final : $t_f = 1.5 T^{100}_{\min}$, où $T^{100}_{\min}$ correspond au temps minimal solution du problème précédent pour $F_\mathrm{max} = 100$.
-**Remarques importantes.**
-* On ne considèrera que des extrémales normales, c'est-à-dire $p^0 = -1$.
-* Le problème a de la structure, le contrôle optimal est composé d'arcs où la poussée est nulle $u(t) = 0$ et où la poussée est maximale $\Vert u(t)\Vert = \gamma_\mathrm{max}$.
-* Vous devrez dans un premier temps résoudre le problème avec le coût suivant ce qui permet de le régulariser et ainsi d'utiliser une méthode de tir simple :
-$$
-\int_{0}^{t_f} \left( \Vert u(t)\Vert - \varepsilon \left( \ln(\Vert u(t)\Vert) + \ln(\gamma_\mathrm{max}-\Vert u(t)\Vert) \right) \right) \mathrm{d}t.
-$$
-%% Cell type:markdown id: tags:
-✏️ **_Question 8:_** Résoudre via du tir simple le problème régularisé pour différentes valeurs de $\varepsilon$ suffisamment petites pour déterminer la structure optimale, en commençant avec $\varepsilon$ suffisamment grand pour faire converger la méthode de tir plus facilement.
-%% Cell type:code id: tags:
-``` julia
-```
-%% Cell type:markdown id: tags:
-✏️ **_Question 9:_** Résoudre via du tir multiple le problème de transfert orbital à consommation minimale. Vous utilisez la structure obtenue précédemment et vous donnerez au solveur une bonne condition initiale à l'aide de la question précédente.
-%% Cell type:code id: tags:
-``` julia
-```
No results found