Skip to content
Snippets Groups Projects
Commit 4498d1af authored by Olivier Cots's avatar Olivier Cots
Browse files

up direct-indirect

parent 32c8dd6c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
tp/direct-indirect.ipynb
+ 18
− 18
View file @ 4498d1af
%% Cell type:markdown id: tags: %% 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) [<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)
# Méthode directe et indirecte # Méthode directe et indirecte
- Date : 2025 - Date : 2025
- Durée approximative : 1.5 séance - Durée approximative : 1.5 séance
Le but est de résoudre par une méthode de tir, 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. Le but est de résoudre par une méthode de tir, 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: %% Cell type:code id: tags:
``` julia ``` julia
using JuMP, Ipopt # pour la méthode directe using JuMP, Ipopt # pour la méthode directe
using DifferentialEquations, NLsolve, ForwardDiff # pour la méthode indirecte using DifferentialEquations, NLsolve, ForwardDiff # pour la méthode indirecte
using Plots # pour les graphiques using Plots # pour les graphiques
include("utils.jl"); # fonctions utilitaires include("utils.jl"); # fonctions utilitaires
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div style="width:95%; <div style="width:95%;
margin:10px; margin:10px;
padding:8px; padding:8px;
color:white; color:white;
background-color: rgb(46, 109, 4); background-color: rgb(46, 109, 4);
border-radius:10px; border-radius:10px;
font-weight:bold; font-weight:bold;
font-size:1.5em; font-size:1.5em;
text-align:center;"> text-align:center;">
Introduction Introduction
</div> </div>
On considère le problème de contrôle optimal suivant : On considère le problème de contrôle optimal suivant :
$$ $$
\left\{ \left\{
\begin{array}{l} \begin{array}{l}
\displaystyle \min_{x, u} \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \\[1.0em] \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] \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} \end{array}
\right. \right.
$$ $$
où $t_f$ = 2. où $t_f$ = 2.
✏️ **Exercice 1.** ✏️ **Exercice 1.**
1. Appliquez le PMP. Que pouvez-vous dire du contrôle maximisant ? 1. Calculer le contrôle maximisant donné par le principe du maximum de Pontryagin (on pourra l'écrire comme une [fonction multivaluée](https://fr.wikipedia.org/wiki/Fonction_multivaluée)) .
2. Peut-on appliquer simplement la méthode de tir simple vu au TP précédent ? 2. Peut-on appliquer simplement la méthode de tir simple vu au TP précédent ?
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div style="width:95%; <div style="width:95%;
margin:10px; margin:10px;
padding:8px; padding:8px;
color:white; color:white;
background-color: rgb(46, 109, 4); background-color: rgb(46, 109, 4);
border-radius:10px; border-radius:10px;
font-weight:bold; font-weight:bold;
font-size:1.5em; font-size:1.5em;
text-align:center;"> text-align:center;">
Méthode directe Méthode directe
</div> </div>
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. Avant de définir la méthode directe, on propose une réécriture du problème. Ce n'est pas une obligation mais cela simplifie l'écriture.
✏️ **Exercice 2.** ✏️ **Exercice 2.**
- Mettez le problème sous forme de Mayer (cf. cours). Vous nommerez la nouvelle variable d'état associée au coût $c(\cdot)$. - Mettre le problème sous forme de Mayer (cf. cours). Vous nommerez la nouvelle variable d'état associée au coût $c(\cdot)$.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Description de la méthode directe.** **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. 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 : 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. 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 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\}\} \{ (x_i, u_i, c_i) ~|~ i \in \{1, \dots, N+1\}\}
$$ $$
l'ensemble des variables d'optimisation du problème discrétisé. 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 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\}. 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é : ✏️ ️**Exercice 3.** Définir pour le problème discrétisé :
1. l'objectif à minimiser en fonction d'une ou plusieurs composantes de $X$ bien choisies. 1. L'objectif à minimiser.
2. les contraintes d'inégalités associées au contrôle. 2. Les contraintes d'inégalités associées au contrôle.
3. les contraintes initiales et finales associées à la variable d'état. 3. Les contraintes initiales et finales associées à la variable d'état et la contrainte de condition initiale pour le coût.
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)). 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: %% Cell type:markdown id: tags:
✏️ ️**Exercice 4.** ✏️ ️**Exercice 4.**
- Modifiez le code suivant afin de résoudre le problème par la méthode directe que vous venez de décrire. - Modifier 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. **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: %% Cell type:code id: tags:
``` julia ``` julia
# Création du modèle JuMP, Utilisation de Ipopt comme solver # Création du modèle JuMP, Utilisation de Ipopt comme solver
sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5)) sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5))
set_optimizer_attribute(sys,"tol",1e-8) set_optimizer_attribute(sys,"tol",1e-8)
set_optimizer_attribute(sys,"constr_viol_tol",1e-6) set_optimizer_attribute(sys,"constr_viol_tol",1e-8)
set_optimizer_attribute(sys,"max_iter",1000) set_optimizer_attribute(sys,"max_iter",1000)
####### #######
####### DEBUT A MODIFIER ####### DEBUT A MODIFIER
####### #######
####### Les ... sont à remplacer ! ####### Les ... sont à remplacer !
####### #######
####### Attention, on écrit le problème sous la forme d'un problème de Mayer ####### Attention, on écrit le problème sous la forme d'un problème de Mayer
####### #######
# Paramètres # Paramètres
t0 = 0 # temps initial t0 = 0 # temps initial
tf = 0 # temps final tf = 0 # temps final
c0 = 0 # coût initial c0 = 0 # coût initial
x0 = 0 # état initial x0 = 0 # état initial
xf = 0 # état final xf = 0 # état final
N = 50 # taille de la grille N = 50 # taille de la grille
Δt = (tf-t0)/N # pas de temps Δt = (tf-t0)/N # pas de temps
@variables(sys, begin @variables(sys, begin
... # coût ... # coût
... # état ... # état
... # control ... # control
end) end)
# Objectif # Objectif
@objective(sys, Min, ...) @objective(sys, Min, ...)
# Conditions initiales et finales # Conditions initiales et finales
@constraints(sys, begin @constraints(sys, begin
con_c0, ... # contraint sur le coût initial con_c0, ... # contraint sur le coût initial
con_x0, ... # contraint sur l'état initial con_x0, ... # contraint sur l'état initial
con_xf, ... # contraint sur l'état final con_xf, ... # contraint sur l'état final
end) end)
# Contraintes de dynamique avec le schéma de Crank-Nicolson # Contraintes de dynamique avec le schéma de Crank-Nicolson
@NLconstraints(sys, begin @NLconstraints(sys, begin
con_dc[j=1:N], ... # contraint différentielle sur le coût con_dc[j=1:N], ... # contraint différentielle sur le coût
con_dx[j=1:N], ... # contraint différentielle sur l'état con_dx[j=1:N], ... # contraint différentielle sur l'état
end); end);
####### #######
####### FIN A MODIFIER ####### FIN A MODIFIER
####### #######
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Solve for the control and state # Solve for the control and state
println("Solving...") println("Solving...")
optimize!(sys) optimize!(sys)
println() println()
# Display results # Display results
if termination_status(sys) == MOI.OPTIMAL if termination_status(sys) == MOI.OPTIMAL
println(" Solution is optimal") println(" Solution is optimal")
elseif termination_status(sys) == MOI.LOCALLY_SOLVED elseif termination_status(sys) == MOI.LOCALLY_SOLVED
println(" (Local) solution found") println(" (Local) solution found")
elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys) 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") println(" Solution is suboptimal due to a time limit, but a primal solution is available")
else else
error(" The model was not solved correctly.") error(" The model was not solved correctly.")
end end
println(" objective value = ", objective_value(sys)) println(" objective value = ", objective_value(sys))
println() println()
# Retrieves values (including duals) # Retrieves values (including duals)
c = value.(c)[:] c = value.(c)[:]
x = value.(x)[:] x = value.(x)[:]
u = value.(u)[:] u = value.(u)[:]
t = (0:N) * value.(Δt) t = (0:N) * value.(Δt)
pc0 = dual(con_c0) pc0 = dual(con_c0)
px0 = dual(con_x0) px0 = dual(con_x0)
pxf = dual(con_xf) pxf = dual(con_xf)
if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end
if(px0*dual(con_dx[1])<0); px0 = -px0; end if(px0*dual(con_dx[1])<0); px0 = -px0; end
if(pxf*dual(con_dx[N])<0); pxf = -pxf; end if(pxf*dual(con_dx[N])<0); pxf = -pxf; end
if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle
sign = -1.0 sign = -1.0
else else
sign = 1.0 sign = 1.0
end end
pc = [ dual(con_dc[i]) for i in 1:N ] pc = [ dual(con_dc[i]) for i in 1:N ]
px = [ dual(con_dx[i]) for i in 1:N ] px = [ dual(con_dx[i]) for i in 1:N ]
pc = sign * [pc0; pc[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 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: %% Cell type:code id: tags:
``` julia ``` julia
x_plot = plot(t, x, xlabel = "t", ylabel = "x", legend = false) x_plot = plot(t, x, xlabel = "t", ylabel = "x", legend = false)
u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre) u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre)
px_plot = plot(t, px, xlabel = "t", ylabel = "p", legend = false) px_plot = plot(t, px, xlabel = "t", ylabel = "p", legend = false)
display(plot(x_plot, px_plot, layout = (1,2), size=(800,300))) display(plot(x_plot, px_plot, layout = (1,2), size=(800,300)))
display(u_plot) display(u_plot)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
✏️ ️**Exercice 5.** ✏️ ️**Exercice 5.**
1. Commentez les résultats. 1. Commenter les résultats.
2. Modifiez la tolérance de l'optimiseur ainsi que le nombre de points de discrétisation. Commentaires. 2. Modifier 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. 3. Décrire la structure du contrôle optimal en fonction du temps.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div style="width:95%; <div style="width:95%;
margin:10px; margin:10px;
padding:8px; padding:8px;
color:white; color:white;
background-color: rgb(46, 109, 4); background-color: rgb(46, 109, 4);
border-radius:10px; border-radius:10px;
font-weight:bold; font-weight:bold;
font-size:1.5em; font-size:1.5em;
text-align:center;"> text-align:center;">
Méthode de tir indirect Méthode de tir indirect
</div> </div>
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. 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. La fonction de commutation est la fonction qui détermine la valeur du contrôle en fonction du signe de celle-ci.
✏️ ️️**Exercice 6.** ✏️ ️️**Exercice 6.**
1. Donner la fonction de commutation. 1. Donner la fonction de commutation.
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'est-à-dire le long d'un arc où la fonction de commutation est constante égale à 0. 2. En dérivant deux fois par rapport au temps la fonction de commutation, 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'est-à-dire le long d'un arc où la fonction de commutation est constante égale à 0.
3. Remplir le code ci-dessous: il faut fournir les 3 contrôles pour définir les flots, puis écrire la fonction de tir. 3. Remplir le code ci-dessous: il faut fournir les 3 contrôles pour définir les flots, puis écrire la fonction de tir.
**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. **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.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# système pseudo-hamiltonien # système pseudo-hamiltonien
function hv(x, p, u) function hv(x, p, u)
return [u, 2x] return [u, 2x]
end end
# systèmes hamiltoniens # systèmes hamiltoniens
hv_min(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur hv_min(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
hv_max(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur hv_max(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
hv_sin(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur hv_sin(x, p) = hv(x, p, NaN) ## REMPLACER NaN par la bonne valeur
# flots # flots
fmin = Flow(hv_min) fmin = Flow(hv_min)
fmax = Flow(hv_max) fmax = Flow(hv_max)
fsin = Flow(hv_sin) fsin = Flow(hv_sin)
# fonction de tir # fonction de tir
function shoot(p0, t1, t2) function shoot(p0, t1, t2)
# integration # integration
x1, p1 = fmin(t0, x0, p0, t1) # x1, p1 sont des scalaires x1, p1 = fmin(t0, x0, p0, t1) # x1, p1 sont des scalaires
x2, p2 = fsin(t1, x1, p1, t2) x2, p2 = fsin(t1, x1, p1, t2)
x3, p3 = fmax(t2, x2, p2, tf) x3, p3 = fmax(t2, x2, p2, tf)
# conditions # conditions
s = zeros(eltype(p0), 3) s = zeros(eltype(p0), 3)
s[1] = NaN ## REMPLACER NaN par la bonne expression s[1] = NaN ## REMPLACER NaN par la bonne expression
s[2] = NaN ## REMPLACER NaN par la bonne expression s[2] = NaN ## REMPLACER NaN par la bonne expression
s[3] = NaN ## REMPLACER NaN par la bonne expression s[3] = NaN ## REMPLACER NaN par la bonne expression
return s return s
end; end;
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
4. Trouver (remplir le code ci-dessous) à partir de la solution de la méthode directe, une bonne initialisation pour la fonction de tir. 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 id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# itéré initiale pour la méthode indirecte de tir # itéré initiale pour la méthode indirecte de tir
p0 = NaN ## REMPLACER NaN par une bonne valeur p0 = NaN ## REMPLACER NaN par une bonne valeur
t1 = NaN ## REMPLACER NaN par une bonne valeur t1 = NaN ## REMPLACER NaN par une bonne valeur
t2 = NaN ## REMPLACER NaN par une bonne valeur t2 = NaN ## REMPLACER NaN par une bonne valeur
# #
y = [ p0 ; t1 ; t2] y = [ p0 ; t1 ; t2]
println("Itéré initial:\n", y) println("Itéré initial:\n", y)
``` ```
%% Output %% Output
Itéré initial: Itéré initial:
[NaN, NaN, NaN] [NaN, NaN, NaN]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# fonction de tir et sa jacobienne # fonction de tir et sa jacobienne
foo(y) = shoot(y...) foo(y) = shoot(y...)
jfoo(y) = ForwardDiff.jacobian(foo, y) jfoo(y) = ForwardDiff.jacobian(foo, y)
# Résolution de shoot(p0, t1, t2) = 0. # Résolution de shoot(p0, t1, t2) = 0.
nl_sol = nlsolve(foo, jfoo, y; xtol=1e-8, method=:trust_region, show_trace=true); nl_sol = nlsolve(foo, jfoo, y; xtol=1e-8, method=:trust_region, show_trace=true);
# Retrieves solution # Retrieves solution
if converged(nl_sol) if converged(nl_sol)
p0 = nl_sol.zero[1] p0 = nl_sol.zero[1]
t1 = nl_sol.zero[2] t1 = nl_sol.zero[2]
t2 = nl_sol.zero[3]; t2 = nl_sol.zero[3];
println("\nFinal solution:\n", nl_sol.zero) println("\nFinal solution:\n", nl_sol.zero)
else else
error("Not converged") error("Not converged")
end end
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# affichage de la solution # affichage de la solution
ode_sol = fmin((t0, t1), x0, p0) ode_sol = fmin((t0, t1), x0, p0)
tt0 = ode_sol.t tt0 = ode_sol.t
xx0 = [ ode_sol[1, j] for j in 1:size(tt0, 1) ] xx0 = [ ode_sol[1, j] for j in 1:size(tt0, 1) ]
pp0 = [ ode_sol[2, j] for j in 1:size(tt0, 1) ] pp0 = [ ode_sol[2, j] for j in 1:size(tt0, 1) ]
uu0 = -ones(size(tt0, 1)) uu0 = -ones(size(tt0, 1))
ode_sol = fsin((t1, t2), xx0[end], pp0[end]) ode_sol = fsin((t1, t2), xx0[end], pp0[end])
tt1 = ode_sol.t tt1 = ode_sol.t
xx1 = [ ode_sol[1, j] for j in 1:size(tt1, 1) ] xx1 = [ ode_sol[1, j] for j in 1:size(tt1, 1) ]
pp1 = [ ode_sol[2, j] for j in 1:size(tt1, 1) ] pp1 = [ ode_sol[2, j] for j in 1:size(tt1, 1) ]
uu1 = zeros(size(tt1, 1)) uu1 = zeros(size(tt1, 1))
ode_sol = fmax((t2, tf), xx1[end], pp1[end]) ode_sol = fmax((t2, tf), xx1[end], pp1[end])
tt2 = ode_sol.t tt2 = ode_sol.t
xx2 = [ ode_sol[1, j] for j in 1:size(tt2, 1) ] xx2 = [ ode_sol[1, j] for j in 1:size(tt2, 1) ]
pp2 = [ ode_sol[2, j] for j in 1:size(tt2, 1) ] pp2 = [ ode_sol[2, j] for j in 1:size(tt2, 1) ]
uu2 = ones(size(tt2, 1)) uu2 = ones(size(tt2, 1))
t_shoot = [ tt0 ; tt1 ; tt2 ] t_shoot = [ tt0 ; tt1 ; tt2 ]
x_shoot = [ xx0 ; xx1 ; xx2 ] x_shoot = [ xx0 ; xx1 ; xx2 ]
p_shoot = [ pp0 ; pp1 ; pp2 ] p_shoot = [ pp0 ; pp1 ; pp2 ]
u_shoot = [ uu0 ; uu1 ; uu2 ] u_shoot = [ uu0 ; uu1 ; uu2 ]
x_plot = plot(t_shoot, x_shoot, xlabel = "t", ylabel = "x", legend = false) x_plot = plot(t_shoot, x_shoot, xlabel = "t", ylabel = "x", legend = false)
u_plot = plot(t_shoot, u_shoot, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre) u_plot = plot(t_shoot, u_shoot, xlabel = "t", ylabel = "u", legend = false, size=(800,300), linetype=:steppre)
px_plot = plot(t_shoot, p_shoot, xlabel = "t", ylabel = "px", legend = false) px_plot = plot(t_shoot, p_shoot, xlabel = "t", ylabel = "px", legend = false)
display(plot(x_plot, px_plot, layout = (1,2), size=(800, 300))) display(plot(x_plot, px_plot, layout = (1,2), size=(800, 300)))
display(u_plot) display(u_plot)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
5. Comparer avec l'affichage de la solution de la méthode directe. 5. Comparer avec l'affichage de la solution de la méthode directe.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment