From 4ff826b513bd88c1322787a27595db60bf318573 Mon Sep 17 00:00:00 2001 From: Oussamamouhtal Date: Thu, 31 Jul 2025 21:29:23 -0400 Subject: [PATCH 1/2] rapport labo 5 --- rapport.qmd | 169 +++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 162 insertions(+), 7 deletions(-) diff --git a/rapport.qmd b/rapport.qmd index b5afcd5..172afba 100644 --- a/rapport.qmd +++ b/rapport.qmd @@ -2,8 +2,8 @@ title: "Rapport de laboratoire 5" subtitle: "MTH8408" author: - - name: Votre nom - email: votre.adresse@polymtl.ca + - name: Mouhtal Oussama + email: oussama-2.mouhtal@polymtl.ca affiliation: - name: Polytechnique Montréal format: @@ -24,45 +24,200 @@ engine: julia --- ```{julia} +#| echo: false #| output: false VERSION == v"1.11.5" || error("please use julia version 1.11.5") using Pkg Pkg.activate("labo9_env") +Pkg.activate("intro-opt-control") +Pkg.add("OptimalControl") +Pkg.add("NLPModelsIpopt") +Pkg.add("Plots") + +using OptimalControl +using NLPModelsIpopt +using Plots ``` # Question 1 -Répondre à l'exercice 1 du laboratoire 8. +Soit le problème suivant : + +\begin{align*} + \min_x \ & \int_0^1 f(x(t), \dot{x}(t), t) \, \mathrm{d}t \\ + \text{s.t.} \ & \int_0^1 h(x(t), \dot{x}(t), t) \, \mathrm{d}t = 0 \\ + & x(0) = x_0, \ x(1) = x_1, +\end{align*} + +Posons $\dot{y}(t) = h(x(t), \dot{x}(t), t)$, la contrainte intégrale devient +\begin{equation*} + \int_0^1 h(x(t), \dot{x}(t), t) \, \mathrm{d}t = \int_0^1 \dot{y}(t)\, \mathrm{d}t = y(1) - y(0) +\end{equation*} +Le problème devient, +\begin{align*} + \min_{x, y} \ & \int_0^1 f(x(t), \dot{x}(t), t) \, \mathrm{d}t \\ + \text{s.t.} \ & \dot{y}(t) = h(x(t), \dot{x}(t), t) \\ + x(0) = x_0, \ & x(1) = x_1, \ y(1) = y(0), +\end{align*} +Pour retrouver un problème de commande optimale, on pose $\dot{x}(t) = u(t)$, +\begin{align*} + \min_{x,y,u} \ & \int_0^1 f(x(t), u(t), t) \, \mathrm{d}t \\ + \text{s.t.} \ \dot{x}(t) = u(t) ,\ & \dot{y}(t) = h(x(t), u(t), t) \\ + x(0) = x_0, \ & x(1) = x_1, \ y(1) = y(0), +\end{align*} + +## Remarque + +La variable d'état $y(t)$ n'agit qu'à partir de sa dérivée dans le problème, donc une translation de $$y(t)$ n'affecte pas sa dérivée, ce qui nous permet de fixer $y(0)$ par exemple à $y(0) =0 $, le problème se simplifie donc à : +\begin{align*} + \min_{x,y,u} \ & \int_0^1 f(x(t), u(t), t) \, \mathrm{d}t \\ + \text{s.t.} \ \dot{x}(t) = u(t) \, & \dot{y}(t) = h(x(t), u(t), t) \\ + x(0) = x_0, \ & x(1) = x_1, \ y(1) = y(0) = 0, +\end{align*} # Question 2 #### Répondre à l'exercice 2 du laboratoire 8 (modélisation et résolution du problème de la corde suspendue). +## Modélisation du problème +L'énergie potentiel de la corde est minimale lorsqu'elle se stabilise. +Le problème peut se modéliser sous cette forme avec contraintes isopérimitrique : + +\begin{align*} + \min_{x} \quad & \rho g \int_0^1 x(t) \sqrt{1 + \dot{x}(t)^2} \, dt \\ +\text{s.t} \quad +& \int_0^1 \sqrt{1 + \dot{x}(t)^2} \, dt = L,\\ +& x(0) = a, \quad x(1) = b. +\end{align*} + +Utilisons le résultat de la première question, on obtient : +\begin{align*} + \min_{x \, y \,u} \quad & \rho g \int_0^1 x(t) \sqrt{1 + u(t)^2} \, dt \\ +\text{s.t} \quad +& \dot{x}(t) = u(t) \, \dot{y}(t) = \sqrt{1 + u(t)^2} - L\\ +& x(0) = a, \quad x(1) = b, \\ +& y(1) = y(0) = 0. +\end{align*} ```{julia} -# votre code ici +t0 = 0 +tf = 1 +g = 9.8 # gravité +x0 = a = 1.0 # hauteur à t=0 +xf = b = 2.0 # hauteur à t=1 +y0 = yf = 0.0 +L = 2.5 # longueur de la corde +ρ = 1.0 # densité + +ocp = @def begin + + t ∈ [ t0, tf ], time + x ∈ R², state + u ∈ R, control + + x(t0) == [x0, y0] + x(tf) == [xf, yf] + + ẋ(t) == [u(t), sqrt(1 + u(t)^2) - L] + + ∫( ρ * g * x1(t)sqrt(1 + u(t)^2) ) → min # énergie potentielle + +end +``` + +```{julia} +sol = solve(ocp) ``` #### Afficher graphiquement les états, commandes et états adjoints finaux. ```{julia} -# votre code ici +plt = plot(sol) +savefig(plt, "solution.png") ``` +![](solution.png) #### Commenter les résultats +- Premièrement, la corde prend une forme parabolique, ce qui paraît intuitif. De même, la fonction $p_1(t)$ est une fonction non constante, comme vu en cours. + +- En ce qui concerne le temps d'exécution, le solveur utilisé a nécessité beaucoup d'itérations pour converger. On observe une convergence superlinéaire lors des dernières itérations. + # Questions 3 #### Modéliser et résoudre numériquement le problème des réservoirs du devoir 5. + +\begin{align*} +\max_{x_1,\, x_2,\, u} \quad & x_2(1) \\ +\text{s.t} \quad +& \dot{x}_1(t) = -x_1(t) + u(t), \\ +& \dot{x}_2(t) = x_1(t), \\ +& x_1(0) = 0, \quad x_2(0) = 0, \\ +& x_1(1) = \frac{1}{2}, \\ +& 0 \leq u(t) \leq 1, \quad \forall t \in [0, 1]. +\end{align*} + ```{julia} -# votre code ici +t0 = 0.0 +tf = 1.0 + +ocp = @def begin + t ∈ [t0, tf], time + x ∈ R², state + u ∈ R, control + + x(t0) == [0.0, 0.0] # L'état des réservoires à t = 0 + ẋ(t) == [-x1(t) + u(t), x1(t) ] + x1(tf) == 0.5 # L'état du réservoires à t = 0 + 0 ≤ u(t) ≤ 1 # Control + + x2(tf) → max +end +``` + + +```{julia} +sol = solve(ocp) ``` #### Afficher graphiquement les états, commandes et états adjoints finaux. ```{julia} -# votre code ici +plt = plot(sol) +savefig(plt, "solution2.png") ``` +![](solution2.png) #### Commenter les résultats +- Le problème suivant est de nature \emph{bang-bang}. Ce que je remarque, c’est que le réservoir 1 dépasse, en un temps inférieur à 1, la moitié du réservoir. À partir de ce moment, le contrôle devient nul afin que le réservoir soit à moitié plein à $t_f = 1$. + +- Dans ce deuxième problème, j’ai imposé que le réservoir ne soit pas rempli à plus de la moitié avant $t_f$. On obtient alors un contrôle non lisse, avec de fortes variations, qui tente de se compenser pour que le réservoir atteigne la moitié de sa capacité à l’instant final. + + + +```{julia} +t0 = 0.0 +tf = 1.0 + +ocp = @def begin + t ∈ [t0, tf], time + x ∈ R², state + u ∈ R, control + + x(t0) == [0.0, 0.0] # L'état des réservoires à t = 0 + ẋ(t) == [-x1(t) + u(t), x1(t) ] + x1(tf) == 0.5 # L'état du réservoires à t = 0 + x1(t) ≤ 0.5 + 0 ≤ u(t) ≤ 1 # Control + + x2(tf) → max +end +sol = solve(ocp) +``` + +```{julia} +plt = plot(sol) +savefig(plt, "solution3.png") +``` +![](solution3.png) From 933db73d4bffcd8257f3e68465c6520394150ca1 Mon Sep 17 00:00:00 2001 From: Dominique Date: Mon, 18 Aug 2025 15:40:30 -0400 Subject: [PATCH 2/2] Update rapport.qmd --- rapport.qmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rapport.qmd b/rapport.qmd index 172afba..e877204 100644 --- a/rapport.qmd +++ b/rapport.qmd @@ -26,7 +26,7 @@ engine: julia ```{julia} #| echo: false #| output: false -VERSION == v"1.11.5" || error("please use julia version 1.11.5") +#VERSION == v"1.11.5" || error("please use julia version 1.11.5") using Pkg Pkg.activate("labo9_env") Pkg.activate("intro-opt-control")