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

foo

parent a646fab2
Branches
No related tags found
No related merge requests found
%% 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)
# Intégration numérique
- Date : 2023
- Durée approximative : 1 séance
%% Cell type:markdown id: tags:
## Introduction
Pour la documentation du package `Julia` sur l'intégration numérique, voir
[documentation de DifferentialEquations.jl](https://diffeq.sciml.ai/stable/).
Il nous faut tout d'abord importer les bons packages de Julia.
%% Cell type:code id: tags:
``` julia
using DifferentialEquations
using ForwardDiff
using LinearAlgebra
using Plots
include("utils.jl"); # plot_traj!, plot_flow!
```
%% Cell type:markdown id: tags:
## Exemple 1
Nous allons ici résoudre le problème à valeur initiale
$$(IVP)\left\{\begin{array}{l}
\dot{x}_1(t)=x_1(t)+x_2(t)+\sin t\\
\dot{x}_2(t)=-x_1(t)+3x_2(t)\\
x_1(0)=-9/25\\
x_2(0)=-4/25,
\end{array}\right.
$$
- Ecrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\dot{x}(t)=f(t,x(t))$.
- Vérifier que la solution de $(IVP)$ est
$$\begin{align*}
x_1(t)&= (-1/25)(13\sin t+9\cos t)\\
x_2(t)&= (-1/25)(3\sin t+4\cos t).
\end{align*}$$
- Coder la fonction $f$ et exécuter le code ci-dessous. Commentaires ?
%% Cell type:code id: tags:
``` julia
# Remarque :
# un problème IVP est une équation différentielle ordinaire
# (EDO en français, ODE en anglais)
# avec en plus une condition initiale.
"""
Second membre de l'IVP
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps
"""
function exemple1(x, λ, t)
xpoint = similar(x)
# A COMPLETER/MODIFIER
xpoint[:] = zeros(size(x))
return xpoint
end
λ = [] # pas de paramètres
t0 = 0
tf = 10
tspan = (t0, tf) # intervalle d'intégration
x0 = [-9/25, -4/25] # condition initiale
prob = ODEProblem(exemple1, x0, tspan, λ) # définition du problème IVP
sol = solve(prob) # résolution du problème IVP
p1 = plot(sol, label = ["x₁(t)" "x₂(t)"], title = "default options") # affichage de la solution
```
%% Cell type:markdown id: tags:
## Le contrôle du pas
Les bons intégrateurs numériques ont une mise à jour automatique du pas. Pour les
[méthodes de Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta)
par exemple, sur un pas $[t_0, t_{1}]$ on calcule par 2 schémas différents deux solutions approchées à l'instant $t_{1}$ : $x_{1}$ et $\hat{x}_{1}$ dont l'erreur locale est respectivement en $O(h^p)$ et $O(h^{p+1})$ (les erreurs globales sont elles en $O(h^{p-1})$ et $O(h^p)$). Ainsi on peut estimer l'erreur locale du schéma le moins précis à l'aide de la différence $x_{1}-\hat{x}_{1}$. On peut alors estimer le pas de façon à avoir la norme de cette différence inférieure à une tolérance donnée `Tol` ou encore à avoir
$$\frac{\|x_{1}-\hat{x}_{1}\|}{\mathrm{Tol}}<1.$$
En pratique on considère des tolérences absolue et relative composante par composante et on définit : $sc_i = \mathrm{Atol}_i + \max(|x_{0i}|,|x_{1i}|)\, \mathrm{Rtol}_i$. On calcule alors le pas de façon à avoir
$$err = \sqrt{\frac{1}{n}\sum_{i=1}^n\left(\frac{x_{1i}-\hat{x}_{1i}}{sc_i}\right)^2}<1.$$
Référence : [Hairer, Nørsett, Wanner (2008) Solving Ordinary Differential Equations I Nonstiff Problems](https://link.springer.com/book/10.1007/978-3-540-78862-1).
Dans `Julia`, on peut modifier les valeurs par défaut : reltol = 1.e-3 et abstol = 1.e-6
(lorsque ces valeurs sont des scalaires, on considère les mêmes quantités pour toutes les composantes).
Réaliser l'intégration numérique pour
- reltol = 1e-3 et abstol = 1e-6
- reltol = 1e-6 et abstol = 1e-9
- reltol = 1e-10 et abstol = 1e-15
et afficher les résultats en ajoutant la solution.
%% Cell type:code id: tags:
``` julia
# calcul avec les options par défaut
sol = solve(prob, reltol = 1e-3, abstol = 1e-6)
p1 = plot(sol, label = ["x₁(t)" "x₂(t)"], title = "reltol=1e-3, abstol=1e-6")
# A COMPLETER/MODIFIER
#
p2 = plot()
#
p3 = plot()
# solution exacte
T = t0:(tf-t0)/100:tf
sinT = sin.(T) # opération vectorielle
cosT = cos.(T)
p4 = plot(T,[(-1/25)*(13*sinT+9*cosT) (-1/25)*(3*sinT+4*cosT)], label = ["x₁(t)" "x₂(t)"], title = "Solution exacte")
# affichage
plot(p1, p2, p3, p4, size=(900, 600))
```
%% Cell type:markdown id: tags:
## Pendule simple
On s'intéresse ici au pendule simple. Les principes physiques de la mécanique classique donnent comme équation qui régit l'évolution du mouvement
$$ ml^2\ddot{\alpha}(t)+mlg\sin(\alpha(t)) +k\dot{\alpha}(t)=0,$$
où $\ddot{\alpha}(t)$ désigne la dérivée seconde de l'angle $\alpha$ par rapport au temps $t$.
- En prenant comme variable d'état $x=(x_1,x_2)=(\alpha, \dot{\alpha})$, écrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\dot{x}(t) = f(t,x(t))$.
- Coder cette fonction ci-dessous et exécuter le code. D'après-vous, observez-vous une oscillation ou une rotation du pendule ?
- Remplacer la vitesse angulaire initiale $\dot{\alpha}(0)$ par $2$. Commentaires.
%% Cell type:code id: tags:
``` julia
"""
Second membre de l'IVP
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.
"""
function pendule(x, λ, t)
xpoint = similar(x)
g, l, k, m = λ
# A COMPLETER/MODIFIER
xpoint[:] = zeros(size(x))
return xpoint
end
#
g = 9.81
l = 10
k = 0
m = 1
λ = [g, l, k, m] # paramètres constants
theta0 = pi/3
x0 = [theta0, 0] # état initial
t0 = 0
tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période
tspan = (t0, tf) # instants initial et terminal
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8) # définition du problème en Julia
sol = solve(prob) # intégration numérique
plot(sol, label = ["x₁(t)" "x₂(t)"]) # affichage de la solution
```
%% Cell type:markdown id: tags:
**Diagramme de phases.**
- Exécutez le code ci-dessous et interprétez le graphique: où sont les oscillations, les rotations et les points d'équilibre stables et instables ?
- On considère le cas où $k=0.15$ (on introduit un frottement). Que se passe-t-il ?
%% Cell type:code id: tags:
``` julia
#
g = 9.81
l = 1.5
k = 0
m = 1
λ = [g, l, k, m] # paramètres constants
plt = plot(xlabel = "x₁", ylabel = "x₂", legend = false) # initialisation du plot
for theta0 in 0:(2*pi)/10:2*pi
theta0_princ = theta0
tf = 3*pi*sqrt(l/g)*(1 + theta0_princ^2/16 + theta0_princ^4/3072) # 2*approximation of the period
tspan = (0.0,tf)
x0 = [theta0 0]
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(plt, sol, vars=(1,2), color="blue") # lw = linewidth
end
theta0 = pi-10*eps()
x0 = [theta0 0]
tf = 50 # problem for tf=50 1/4 of the period!
tspan = (0.0,tf)
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(plt, sol, vars=(1,2), xlims = (-2*pi,4*pi), color="green") # lw = linewidth
theta0 = pi+10*eps()
x0 = [theta0 0]
tf = 50
tspan = (0.0,tf)
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(plt, sol, vars=(1,2), xlims = (-2*pi,4*pi), color="green") # lw = linewidth
# circulation case
for thetapoint0 in 0:1.:4
tf = 10
tspan = (0.,tf)
x0 = [-pi thetapoint0] # thetapoint0 > 0 so theta increases from -pi to ...
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(plt, sol, vars=(1,2), color="red") # lw = linewidth
end
for thetapoint0 in -4:1.:0
tf = 10
tspan = (0.,tf)
x0 = [3*pi thetapoint0] # thetapoint0 < 0 so theta decreases from 3pi to ...
prob = ODEProblem(pendule, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(plt, sol, vars=(1,2), color="purple") # lw = linewidth
end
plot!(plt, [-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter)
plot!(plt, xlims = (-pi,3*pi), size=(900, 600))
```
%% Cell type:markdown id: tags:
## Modèle de Van der Pol
L'équation différentielle considérée est l'équation de Van der Pol
$$(IVP)\left\{\begin{array}{l}
\dot{x}_1(t)=x_2(t)\\
\dot{x}_2(t)=(1-x_1^2(t))x_2(t)-x_1(t)\\
x_1(0)=2.00861986087484313650940188\\
x_2(0)=0
\end{array}\right.
$$
jusqu'au temps $t_f=T=6.6632868593231301896996820305$, où $T$ est la période de la solution.
Codez le second membre de l'IVP ci-dessous et exécutez le code. Commentaires.
%% Cell type:code id: tags:
``` julia
"""
Second membre de l'IVP: modèle de Van der Pol
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.
"""
function vdp(x, λ, t)
xpoint = similar(x)
# A COMPLETER/MODIFIER
xpoint[:] = zeros(size(x))
return xpoint
end
#
t0 = 0
tf = 6.6632868593231301896996820305
tspan = (t0, tf)
#
plot(xlabel = "x₁", ylabel = "x₂", legend = false)
# Trajectoires intérieures
for x01 in -2:0.4:0
x0 = [x01,0]
prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(sol, vars=(1,2), color = "blue")
end
# Trajectoires extérieures
for x02 in 2.5:0.5:4
x0 = [0,x02]
prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(sol, vars=(1,2), color = "green")
end
# Cycle limite périodique
x0 = [2.00861986087484313650940188,0]
prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)
sol = solve(prob)
plot!(sol, vars=(1,2), color = "red", lw = 2.)
#
plot!([0], [0], seriestype=:scatter) # point d'équilibre
plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(900, 600))
```
%% Cell type:markdown id: tags:
**Stabilité.**
On calcule ci-dessous les valeurs propres de la matrice jacobienne du second membre de l'IVP au point d'équilibre $x_e = (0, 0)$. Commentaires (rappelez-vous votre cours d'[automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf)).
%% Cell type:code id: tags:
``` julia
# Jacobienne de la fonction vdp
dvdp(x) = ForwardDiff.jacobian(x -> vdp(x, [], 0), x)
# Point d'équilibre et la matrice Jacobienne en ce point
xe = [0, 0]
A = dvdp(xe)
# Valeurs propres de la matrice Jacobienne
eigvals(A)
```
%% Cell type:markdown id: tags:
## Modèle proie-prédateur (modèle de Lotka-Volterra)
L'équation différentielle considérée est
$$\left\{\begin{array}{l}
\dot{x}_1(t)= \phantom{-}x_1(t) ( \alpha - \beta x_2(t)) \\[0.5em]
\dot{x}_2(t)= -x_2(t) ( \gamma - \delta x_1(t))
\end{array}\right.
$$
- $t$ est le temps ;
- $x(t)$ est l'effectif des proies en fonction du temps ;
- $y(t)$ est l'effectif des prédateurs en fonction du temps ;
Les paramètres suivants caractérisent les interactions entre les deux espèces :
- $\alpha$, taux de reproduction intrinsèque des proies (constant, indépendant du nombre de prédateurs) ;
- $\beta$, taux de mortalité des proies dû aux prédateurs rencontrés ;
- $\delta$, taux de reproduction des prédateurs en fonction des proies rencontrées et mangées ;
- $\gamma$, taux de mortalité intrinsèque des prédateurs (constant, indépendant du nombre de proies) ;
**Questions.**
- Le point $(0, 0)$ est clairement un point d'équilibre stable. Il existe un autre point d'équilibre, lequel ?
- Afficher l'autre point d'équilibre ainsi que quelques trajectoires du système dans le plan de phase, qui rentrent dans les limites du plot ci-dessous. Qu'observez-vous sur les trajectoires.
%% Cell type:code id: tags:
``` julia
# ------------------------------
# DEBUT A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP
# ...
# FIN A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP
# ------------------------------
# paramètres
α = 2/3
β = 4/3
γ = 1
δ = 1
λ = [α, β, γ, δ]
t0 = 0
tf = 20
tspan = (t0, tf)
plt = plot(xlabel = "x₁", ylabel = "x₂", legend = false)
# ------------------------------
# DEBUT A COMPLETER/MODIFIER L'AFFICHAGE
# ...
# FIN A COMPLETER/MODIFIER L'AFFICHAGE
# ------------------------------
plot!(xlims = (0, 4), ylims = (0, 2.5), size=(900, 600))
```
%% Cell type:markdown id: tags:
## Flots en dimension 2
On considère le problème de contrôle optimal suivant:
$$ \frac{1}{2}\int_0^{1} u^2(t)\,\mathrm{d}t \to \min, $$
sous les contraintes
$$ \dot{x}(t) = -x(t) + u(t)$$
$$ x(0)=x_0=-1,\quad x(1)=0. $$
Le but est ici de tracer les extrémales et le flot associé pour des points initiaux $z(0) = (x_0, p(0))$.
%% Cell type:markdown id: tags:
On note $z = (x, p)$ la paire état, état adoint. On note $F(t, z)$ le système hamiltonien donné par le PMP, c-a-d:
$$
F(t, z) = \left(\frac{\partial H}{\partial p}(t, z, u(z)), -\frac{\partial H}{\partial x}(t, z, u(z))\right)
$$
où $u(z)$ est le contrôle maximisant.
- Codez ci-dessous la fonction $F$ et visualiser son flot.
- Visualisez l'extrémale la plus proche de la solution.
- Qu'observez-vous sur la forme des courbes vertes, c-a-d des fronts d'ondes ? Pourquoi ?
**Remarque.** On rappelle que les extrémales s'écrivent $x(t) = p_0\, \mathrm{sh}(t) + x_0\, e^{-t}$, $p(t) = p_0\, e^t$.
%% Cell type:code id: tags:
``` julia
# système hamiltonien
function F(z)
zpoint = similar(z)
# A COMPLETER/MODIFIER
zpoint = zeros(size(z))
return zpoint
end
t0 = 0
tf = 1
# initialisation du plot
plt = plot()
ymin = -0.2
ymax = 2.5
# plot de trajectoires
z0 = [[-1, p0] for p0 range(ymin, ymax, length=50)]
plot_traj!(plt, F, z0, tf)
# plot du flot
z0 = [[-1, p0] for p0 range(ymin, ymax, length=50)]
plot_flow!(plt, F, z0, tf)
# affichages additionels
plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax))
plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)
plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(900, 600), xlabel = "x", ylabel = "p")
```
%% Cell type:markdown id: tags:
On maintenant considère le problème de contrôle optimal suivant:
$$ \frac{1}{2}\int_0^{1} u^2(t)\,\mathrm{d}t \to \min, $$
sous les contraintes
$$ \dot{x}(t) = -x(t) + \alpha x^2(t) + u(t)$$
$$ x(0)=x_0=-1,\quad x(1)=0. $$
Le but est ici de tracer les extrémales et le flot associé pour des points initiaux $z(0) = (x_0, p(0))$. Visualisez le flot et le tp est terminé.
%% Cell type:code id: tags:
``` julia
# paramètre
α = 1.5
# système hamiltonien
function F(z)
zpoint = similar(z)
zpoint[1] = -z[1] + α*z[1]^2 + z[2]
zpoint[2] = (1 - 2*α*z[1])*z[2]
return zpoint
end
t0 = 0
tf = 1
# initialisation du plot
plt = plot()
ymin = -0.2
ymax = 2.5
# plot de trajectoires
z0 = [[-1, p0] for p0 range(ymin, ymax, length=50)]
plot_traj!(plt, F, z0, tf)
# plot du flot
z0 = [[-1, p0] for p0 range(ymin, ymax, length=100)]
plot_flow!(plt, F, z0, tf)
# affichages additionels
plot!(plt, xlims = (-1.2, 2), ylims = (ymin, ymax))
plot!(plt, [-1, -1], [ymin, ymax], linestyle = :dot, color = :black)
plot!(plt, [0, 0], [ymin, ymax], linestyle = :dot, color = :black, size=(900, 600), xlabel = "x", ylabel = "p")
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment