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

foo

parent 561f5b04
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 150 additions and 19 deletions
tp/animation-etu.ipynb
+ 150
− 19
View file @ 31e4e23c
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
![n7](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png) ![n7](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png)
# Transfert orbital en temps minimal # Transfert orbital en temps minimal
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. 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.
---- ----
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 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 \left\lbrace
\begin{array}{l} \begin{array}{l}
\min J(x, u, t_f) = t_f \\[1.0em] \min J(x, u, t_f) = t_f \\[1.0em]
\ \ \dot{x}_{1}(t) = ~ x_{3}(t) \\ \ \ \dot{x}_{1}(t) = ~ x_{3}(t) \\
\ \ \dot{x}_{2}(t) = ~ x_{4}(t) \\ \ \ \dot{x}_{2}(t) = ~ x_{4}(t) \\
\ \ \dot{x}_{3}(t) = -\dfrac{\mu\, x_{1}(t)}{r^{3}(t)} + u_{1}(t) \\ \ \ \dot{x}_{3}(t) = -\dfrac{\mu\, x_{1}(t)}{r^{3}(t)} + u_{1}(t) \\
\ \ \dot{x}_{4}(t) = -\dfrac{\mu\, x_{2}(t)}{r^{3}(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.0em] \ \ \dot{x}_{4}(t) = -\dfrac{\mu\, x_{2}(t)}{r^{3}(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.0em]
\ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\ \ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\
\ \ r^2(t_f) = r_{f}^2, ~~ 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), \\ \ \ r^2(t_f) = r_{f}^2, ~~ 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} \end{array}
\right. \right.
$$ $$
avec $r(t)=\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$. avec $r(t)=\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$.
Les unités choisies sont le kilomètre pour les distances et l'heure pour les temps. On donne les paramètres suivants : 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}. \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 : 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} \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}$. où m est la masse du satellite qu'on fixe à $m=2000\ \mathrm{kg}$.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Résolution via du tir simple indirect ## Résolution via du tir simple indirect
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Packages utiles # Packages utiles
using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Les constantes du pb # Les constantes du pb
x0 = [-42272.67, 0, 0, -5796.72] # état initial x0 = [-42272.67, 0, 0, -5796.72] # état initial
μ = 5.1658620912*1e12 μ = 5.1658620912*1e12
rf = 42165.0 ; rf = 42165.0 ;
F_max = 100.0 F_max = 100.0
γ_max = F_max*3600.0^2/(2000.0*10^3) γ_max = F_max*3600.0^2/(2000.0*10^3)
t0 = 0.0 t0 = 0.0
rf3 = rf^3 ; rf3 = rf^3 ;
α = sqrt(μ/rf3); α = sqrt(μ/rf3);
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Function to get the flow of a Hamiltonian system # Function to get the flow of a Hamiltonian system
function Flow(hv) function Flow(hv)
function rhs!(dz, z, dummy, t) function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2 n = size(z, 1)÷2
dz[:] = hv(z[1:n], z[n+1:2*n]) dz[:] = hv(z[1:n], z[n+1:2*n])
end end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.1) function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.05)
z0 = [ x0 ; p0 ] z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan) ode = ODEProblem(rhs!, z0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat) sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol return sol
end end
function f(t0, x0, p0, t; abstol=1e-12, reltol=1e-12, saveat=[]) function f(t0, x0, p0, t; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, t), x0, p0; abstol=abstol, reltol=reltol, saveat=saveat) sol = f((t0, t), x0, p0; abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1) n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end] return sol[1:n, end], sol[n+1:2*n, end]
end end
return f return f
end; end;
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-info"> <div class="alert alert-info">
**_Question 1:_** **_Question 1:_**
* Donner le pseudo-hamiltonien associé au problème de contrôle optimal. * Donner le pseudo-hamiltonien associé au problème de contrôle optimal.
* Calculer le contrôle maximisant dans le cas normal ($p^0=-1$). On supposera que $(p_3, p_4)\neq (0,0)$). * Calculer le contrôle maximisant dans le cas normal ($p^0=-1$). On supposera que $(p_3, p_4)\neq (0,0)$).
* Donner le pseudo système hamiltonien * Donner le pseudo système hamiltonien
$$ $$
\vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u), \vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-\frac{\partial H}{\partial x}(x, p, u) \right). -\frac{\partial H}{\partial x}(x, p, u) \right).
$$ $$
</div> </div>
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
**Réponse** (double cliquer pour modifier le texte) **Réponse** (double cliquer pour modifier le texte)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
<div class="alert alert-info"> <div class="alert alert-info">
**_Question 2:_** Compléter le code suivant pour calculer la trajectoire optimale et l'afficher. **_Question 2:_** Compléter le code suivant pour calculer la trajectoire optimale et l'afficher.
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
##### #####
##### A COMPLETER ##### A COMPLETER
# Contrôle maximisant # Contrôle maximisant
function control(p) function control(p)
u = zeros(eltype(p),2) u = zeros(eltype(p),2)
u[1] = ... u[1] = ...
u[2] = ... u[2] = ...
return u return u
end; end;
# Hamiltonien maximisé # Hamiltonien maximisé
function hfun(x, p) function hfun(x, p)
h = ... h = ...
return h return h
end end
# Système hamiltonien # Système hamiltonien
function hv(x, p) function hv(x, p)
n = size(x, 1) n = size(x, 1)
hv = zeros(eltype(x), 2*n) hv = zeros(eltype(x), 2*n)
... ...
return hv return hv
end end
f = Flow(hv); f = Flow(hv);
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
On note On note
$ $
\alpha := \sqrt{\frac{\mu}{r_f^3}}. \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$ La condition terminale peut se mettre sous la forme $c(x(t_f)) = 0$, avec $c \colon \mathbb{R}^4 \to \mathbb{R}^3$
donné par donné par
$$ $$
c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \alpha\, x_{2}, ~~ x_{4} - \alpha\, x_{1}). c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \alpha\, x_{2}, ~~ x_{4} - \alpha\, x_{1}).
$$ $$
Le temps final étant libre, on a la condition au temps final 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{(cas normal)} H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(cas normal)}
$$ $$
De plus, la condition de transversalité De plus, la condition de transversalité
$$ $$
p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3, p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
$$ $$
conduit à conduit à
$$ $$
\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. \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 les conditions aux limites, la condition finale sur le pseudo-hamiltonien et la condition de transversalité, la fonction de tir simple est donnée par En considérant les conditions 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{equation*}
\begin{array}{rlll} \begin{array}{rlll}
S \colon & \mathbb{R}^5 & \longrightarrow & \mathbb{R}^5 \\ S \colon & \mathbb{R}^5 & \longrightarrow & \mathbb{R}^5 \\
& (p_0, t_f) & \longmapsto & & (p_0, t_f) & \longmapsto &
S(p_0, t_f) := \begin{pmatrix} S(p_0, t_f) := \begin{pmatrix}
c(x(t_f, x_0, p_0)) \\[0.5em] c(x(t_f, x_0, p_0)) \\[0.5em]
\Phi(z(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 H(z(t_f, x_0, p_0), u(z(t_f, x_0, p_0))) - 1
\end{pmatrix} \end{pmatrix}
\end{array} \end{array}
\end{equation*} \end{equation*}
où $z(t_f, x_0, p_0)$ est la solution au temps de $t_f$ du pseudo système hamiltonien bouclé par où $z(t_f, x_0, p_0)$ est la solution au temps de $t_f$ du pseudo système hamiltonien bouclé par
le contrôle maximisant, partant au temps $t_0$ du point $(x_0, p_0)$. On rappelle que l'on note le contrôle maximisant, partant au temps $t_0$ du point $(x_0, p_0)$. On rappelle que l'on note
$z=(x, p)$. $z=(x, p)$.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
##### #####
##### A COMPLETER ##### A COMPLETER
# Fonction de tir # Fonction de tir
function shoot(p0, tf) function shoot(p0, tf)
s = zeros(eltype(p0), 5) s = zeros(eltype(p0), 5)
... ...
return s return s
end; end;
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Itéré initial # Itéré initial
y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4] # pour F_max = 100N 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 # Jacobienne de la fonction de tir
foo(y) = shoot(y[1:4], y[5]) foo(y) = shoot(y[1:4], y[5])
jfoo(y) = ForwardDiff.jacobian(foo, y) jfoo(y) = ForwardDiff.jacobian(foo, y)
# Résolution de shoot(p0, tf) = 0 # Résolution de shoot(p0, tf) = 0
nl_sol = nlsolve(foo, jfoo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true); nl_sol = nlsolve(foo, jfoo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true);
# On récupère la solution si convergence # On récupère la solution si convergence
if converged(nl_sol) if converged(nl_sol)
p0_sol_100 = nl_sol.zero[1:4] p0_sol_100 = nl_sol.zero[1:4]
tf_sol_100 = nl_sol.zero[5] tf_sol_100 = nl_sol.zero[5]
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
# Corps du satellite
@userplot SatBody
@recipe function f(cp::SatBody)
x, y = cp.args
seriestype --> :shape
fillcolor --> :goldenrod1
linecolor --> :goldenrod1
x, y
end
@userplot SatLine
@recipe function f(cp::SatLine)
x, y = cp.args
linecolor --> :black
x, y
end
# Bras du satellite
@userplot SatBras
@recipe function f(cp::SatBras)
x, y = cp.args
seriestype --> :shape
fillcolor --> :goldenrod1
linecolor --> :white
x, y
end
# panneau
@userplot PanBody
@recipe function f(cp::PanBody)
x, y = cp.args
seriestype --> :shape
fillcolor --> :dodgerblue4
linecolor --> :black
x, y
end
function satellite!(pl; position=[0;0], scale=1, rotate=0)
# Fonctions utiles
R(θ) = [ cos(θ) -sin(θ)
sin(θ) cos(θ)] # rotation
T(x, v) = x.+v # translation
H(λ, x) = λ.*x # homotéthie
SA(x, θ, c) = c .+ 2*[cos(θ);sin(θ)]'*(x.-c).*[cos(θ);sin(θ)]-(x.-c) # symétrie axiale
#
O = position
Rθ = R(rotate)
# Paramètres
α = π/10.0 # angle du tube, par rapport à l'horizontal
β = π/2-π/40 # angle des bras, par rapport à l'horizontal
Lp = scale.*1.4*cos(α) # longueur d'un panneau
lp = scale.*2.6*sin(α) # largueur d'un panneau
# Corps du satellite
t = range(-α, α, length = 50)
x = scale.*cos.(t);
y = scale.*sin.(t);
Δx = scale.*cos(α)
Δy = scale.*sin(α)
M = T(Rθ*[ x'; y'], O); satbody!(pl, M[1,:], M[2,:]) # bord droit
M = T(Rθ*[-x'; y'], O); satbody!(pl, M[1,:], M[2,:]) # bord gauche
M = T(Rθ*[scale.*[-cos(α), cos(α), cos(α), -cos(α)]'
scale.*[-sin(α), -sin(α), sin(α), sin(α)]'], O);
satbody!(pl, M[1,:], M[2,:]) # interieur
M = T(Rθ*[ x'; y'], O); satline!(pl, M[1,:], M[2,:]) # bord droit
M = T(Rθ*[-x'; y'], O); satline!(pl, M[1,:], M[2,:]) # bord gauche
M = T(Rθ*[ x'.-2*Δx; y'], O); satline!(pl, M[1,:], M[2,:]) # bord gauche (droite)
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:]) # haut
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
-scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:]) # bas
# Bras du satellite
lb = scale.*3*cos(α) # longueur des bras
eb = scale.*cos(α)/30 # demi largeur des bras
xb = 0.0
yb = scale.*sin(α)
M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
[yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb]'], O);
satbras!(pl, M[1,:], M[2,:])
M = T(Rθ*[ [xb-eb, xb+eb, xb+eb+lb*cos(β), xb-eb+lb*cos(β), xb-eb]';
-[yb, yb, yb+lb*sin(β), yb+lb*sin(β), yb']'], O);
satbras!(pl, M[1,:], M[2,:])
# panneau
ep = (lb-3*lp)/6
panneau = [0 Lp Lp 0
0 0 lp lp]
ey = 3*eb # eloignement des panneaux au bras
vy = [cos(β-π/2); sin(β-π/2)] .* ey
v0 = [0; yb]
v1 = 2*ep*[cos(β); sin(β)]
v2 = (3*ep+lp)*[cos(β); sin(β)]
v3 = (4*ep+2*lp)*[cos(β); sin(β)]
pa1 = T(R(β-π/2)*panneau, v0+v1+vy); pa = T(Rθ*pa1, O); panbody!(pl, pa[1,:], pa[2,:])
pa2 = T(R(β-π/2)*panneau, v0+v2+vy); pa = T(Rθ*pa2, O); panbody!(pl, pa[1,:], pa[2,:])
pa3 = T(R(β-π/2)*panneau, v0+v3+vy); pa = T(Rθ*pa3, O); panbody!(pl, pa[1,:], pa[2,:])
pa4 = SA(pa1, β, [xb; yb]); pa = T(Rθ*pa4, O); panbody!(pl, pa[1,:], pa[2,:])
pa5 = SA(pa2, β, [xb; yb]); pa = T(Rθ*pa5, O); panbody!(pl, pa[1,:], pa[2,:])
pa6 = SA(pa3, β, [xb; yb]); pa = T(Rθ*pa6, O); panbody!(pl, pa[1,:], pa[2,:])
pa7 = SA(pa1, 0, [0; 0]); pa = T(Rθ*pa7, O); panbody!(pl, pa[1,:], pa[2,:])
pa8 = SA(pa2, 0, [0; 0]); pa = T(Rθ*pa8, O); panbody!(pl, pa[1,:], pa[2,:])
pa9 = SA(pa3, 0, [0; 0]); pa = T(Rθ*pa9, O); panbody!(pl, pa[1,:], pa[2,:])
pa10 = SA(pa7, -β, [xb; -yb]); pa = T(Rθ*pa10, O); panbody!(pl, pa[1,:], pa[2,:])
pa11 = SA(pa8, -β, [xb; -yb]); pa = T(Rθ*pa11, O); panbody!(pl, pa[1,:], pa[2,:])
pa12 = SA(pa9, -β, [xb; -yb]); pa = T(Rθ*pa12, O); panbody!(pl, pa[1,:], pa[2,:])
#
#display(pl)
end;
@userplot TrajectoryPlot @userplot TrajectoryPlot
@recipe function f(cp::TrajectoryPlot) @recipe function f(cp::TrajectoryPlot)
t, x, y, tf = cp.args t, x, y, tf = cp.args
n = argmin(abs.(t.-tf)) n = argmin(abs.(t.-tf))
inds = 1:n inds = 1:n
seriescolor --> :white seriescolor --> :white
linewidth --> range(0, 5, length = n) linewidth --> range(0, 5, length = n)
seriesalpha --> range(0, 1, length = n) seriesalpha --> range(0, 1, length = n)
aspect_ratio --> 1 aspect_ratio --> 1
label --> false label --> false
x[inds], y[inds] x[inds], y[inds]
end end
# Fonction d'affichage d'une solution # Fonction d'affichage d'une solution
function plot_solution(p0, tf; animation=false, fps=15, nFrame=100) function plot_solution(p0, tf; animation=false, fps=15, nFrame=200)
# On trace l'orbite de départ et d'arrivée # On trace l'orbite de départ et d'arrivée
gr(dpi=300, size=(500,400), thickness_scaling=1) gr(dpi=300, size=(500,400), thickness_scaling=1)
r0 = norm(x0[1:2]) r0 = norm(x0[1:2])
v0 = norm(x0[3:4]) v0 = norm(x0[3:4])
a = 1.0/(2.0/r0-v0*v0/μ) a = 1.0/(2.0/r0-v0*v0/μ)
t1 = r0*v0*v0/μ - 1.0; t1 = r0*v0*v0/μ - 1.0;
t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ); t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ);
e_ellipse = norm([t1 t2]) e_ellipse = norm([t1 t2])
p_orb = a*(1-e_ellipse^2); p_orb = a*(1-e_ellipse^2);
n_theta = 101 n_theta = 151
Theta = range(0, stop=2*pi, length=n_theta) Theta = range(0, stop=2*pi, length=n_theta)
X1_orb_init = zeros(n_theta) X1_orb_init = zeros(n_theta)
X2_orb_init = zeros(n_theta) X2_orb_init = zeros(n_theta)
X1_orb_arr = zeros(n_theta) X1_orb_arr = zeros(n_theta)
X2_orb_arr = zeros(n_theta) X2_orb_arr = zeros(n_theta)
# Orbite initiale # Orbite initiale
for i in 1:n_theta for i in 1:n_theta
theta = Theta[i] theta = Theta[i]
r_orb = p_orb/(1+e_ellipse*cos(theta)); r_orb = p_orb/(1+e_ellipse*cos(theta));
X1_orb_init[i] = r_orb*cos(theta); X1_orb_init[i] = r_orb*cos(theta);
X2_orb_init[i] = r_orb*sin(theta); X2_orb_init[i] = r_orb*sin(theta);
end end
# Orbite d'arrivée # Orbite d'arrivée
for i in 1:n_theta for i in 1:n_theta
theta = Theta[i] theta = Theta[i]
X1_orb_arr[i] = rf*cos(theta) ; X1_orb_arr[i] = rf*cos(theta) ;
X2_orb_arr[i] = rf*sin(theta); X2_orb_arr[i] = rf*sin(theta);
end; end;
# Calcul de la trajectoire # Calcul de la trajectoire
ode_sol = f((t0, tf), x0, p0) ode_sol = f((t0, tf), x0, p0)
t = ode_sol.t t = ode_sol.t
n = size(t, 1) n = size(t, 1)
x1 = [ode_sol[1, j] for j in 1:n ] x1 = [ode_sol[1, j] for j in 1:n ]
x2 = [ode_sol[2, j] for j in 1:n ] x2 = [ode_sol[2, j] for j in 1:n ]
v1 = [ode_sol[3, j] for j in 1:n ]
v2 = [ode_sol[4, j] for j in 1:n ]
u = zeros(2, length(t)) u = zeros(2, length(t))
for j in 1:size(t, 1) for j in 1:size(t, 1)
u[:,j] = control(ode_sol[5:8, j]) u[:,j] = control(ode_sol[5:8, j])
end end
ee = 0.1 ee = 0.1
xmin = minimum([x1; X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin); xmin = minimum([x1; X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);
xmax = maximum([x1; X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax); xmax = maximum([x1; X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);
ymin = minimum([x2; X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin); ymin = minimum([x2; X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);
ymax = maximum([x2; X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax); ymax = maximum([x2; X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);
if animation if animation
nFrame = min(nFrame, n); nFrame = min(nFrame, n);
anim = @animate for i ∈ 1:nFrame anim = @animate for i ∈ 1:nFrame
px = plot(background_color=:black, xlims=(xmin, xmax), ylims=(ymin, ymax), px = plot(background_color=:gray26, xlims=(xmin, xmax), ylims=(ymin, ymax),
legend = false, border=:none, aspect_ratio=:equal) legend = false, framestyle = :none, aspect_ratio=:equal)
plot!(px, X1_orb_init, X2_orb_init, color = :green, linewidth=2) plot!(px, X1_orb_init, X2_orb_init, color = :olivedrab1, linewidth=2)
plot!(px, X1_orb_arr, X2_orb_arr, color = :blue, linewidth=2) plot!(px, X1_orb_arr, X2_orb_arr, color = :turquoise1, linewidth=2)
plot!(px, [0.0], [0.0], color = :blue, seriestype=:scatter, markersize = 20, markerstrokewidth=0) plot!(px, [0.0], [0.0], color = :steelblue2, seriestype=:scatter, markersize = 20, markerstrokewidth=0)
trajectoryplot!(px, t, x1, x2, i*tf/nFrame) trajectoryplot!(px, t, x1, x2, i*tf/nFrame)
tcur = i*tf/nFrame tcur = i*tf/nFrame
indf = argmin(abs.(t.-tcur)) indf = argmin(abs.(t.-tcur))
plot!(px, [x1[indf]], [x2[indf]], seriestype=:scatter, color = :red, markerstrokewidth=0) plot!(px, [x1[1]], [x2[1]], seriestype=:scatter, color = :white, markerstrokewidth=0)
satellite!(px, position=[x1[indf];x2[indf]], scale=2000, rotate=atan(v2[indf], v1[indf]))
end end
gif(anim, "anim.gif", fps=fps); gif(anim, "anim.gif", fps=fps);
else else
pu1 = plot(t, u[1,:], color = :red, xlabel = "t", ylabel = "u1", legend = false) pu1 = plot(t, u[1,:], color = :red, xlabel = "t", ylabel = "u1", legend = false)
pu2 = plot(t, u[2,:], color = :red, xlabel = "t", ylabel = "u2", legend = false) pu2 = plot(t, u[2,:], color = :red, xlabel = "t", ylabel = "u2", legend = false)
px = plot(background_color=:black, xlims=(xmin, xmax), ylims=(ymin, ymax), px = plot(background_color=:gray26, xlims=(xmin, xmax), ylims=(ymin, ymax),
legend = false, border=:none, aspect_ratio=:equal) legend = false, framestyle = :none, aspect_ratio=:equal)
plot!(px, X1_orb_init, X2_orb_init, color = :green, linewidth=2) plot!(px, X1_orb_init, X2_orb_init, color = :olivedrab1, linewidth=2)
plot!(px, X1_orb_arr, X2_orb_arr, color = :blue, linewidth=2) plot!(px, X1_orb_arr, X2_orb_arr, color = :turquoise1, linewidth=2)
plot!(px, [0.0], [0.0], color = :blue, seriestype=:scatter, markersize = 20, markerstrokewidth=0) plot!(px, [0.0], [0.0], color = :steelblue2, seriestype=:scatter, markersize = 20, markerstrokewidth=0)
plot!(px, [x1[1]], [x2[1]], seriestype=:scatter, color = :red, markerstrokewidth=0) plot!(px, [x1[1]], [x2[1]], seriestype=:scatter, color = :white, markerstrokewidth=0)
plot!(x1, x2, color = :white, size = (800, 800), linewidth = range(0, 5, length = n)) plot!(px, x1, x2, color = :white, size = (800, 800), linewidth = range(0, 5, length = n))
satellite!(px, position=[x1[end];x2[end]], scale=2000, rotate=atan(v2[end], v1[end]))
display(plot(pu1, pu2, layout = (1,2), size = (800,400))) display(plot(pu1, pu2, layout = (1,2), size = (800,400)))
display(px) display(px)
end end
end; end;
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# On affiche la solution pour Fmax = 100 # On affiche la solution pour Fmax = 100
plot_solution(p0_sol_100, tf_sol_100) plot_solution(p0_sol_100, tf_sol_100)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Etude du temps de transfert en fonction de la poussée maximale # Etude du temps de transfert en fonction de la poussée maximale
<div class="alert alert-info"> <div class="alert alert-info">
**_Question 3:_** **_Question 3:_**
* A l'aide de ce que vous avez fait précédemment, déterminer $t_f$ pour * A l'aide de ce que vous avez fait précédemment, déterminer $t_f$ pour
$$ $$
F_\mathrm{max} \in \{100, 90, 80, 70, 60, 50, 40, 30, 20 \}. 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. * Tracer ensuite $t_f$ en fonction de $F_\mathrm{max}$ et commenter le résultat.
</div> </div>
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Les différentes valeurs de poussées # Les différentes valeurs de poussées
F_max_span = range(100, stop=20, length=11); 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)]; γ_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: %% Cell type:code id: tags:
``` julia ``` julia
##### #####
##### A COMPLETER ##### A COMPLETER
# Solution calculée précédemment # Solution calculée précédemment
y_guess = [p0_sol_100; [tf_sol_100]] y_guess = [p0_sol_100; [tf_sol_100]]
tf_sols = zeros(length(γ_max_span)) tf_sols = zeros(length(γ_max_span))
p0_sols = zeros(4, length(γ_max_span)) p0_sols = zeros(4, length(γ_max_span))
... ...
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Affichage de tf en fonction de la poussée maximale # Affichage de tf en fonction de la poussée maximale
plot(F_max_span, tf_sols, aspect_ratio=:equal, legend=false, xlabel="Fmax", ylabel="tf") plot(F_max_span, tf_sols, aspect_ratio=:equal, legend=false, xlabel="Fmax", ylabel="tf")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` julia ``` julia
# Plots sol pour F_max = 20N # Plots sol pour F_max = 20N
plot_solution(p0_sols[:, end], tf_sols[end], animation=false, fps=10) plot_solution(p0_sols[:, end], tf_sols[end], animation=false)
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment