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

foo

parent 383c9000
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
transfert-temps-min-original.mp4
+ 0
− 0
View file @ 0ec9fcf2
No preview for this file type
transfert-temps-min.gif
+ 0
− 0
View replaced file @ 383c9000 View file @ 0ec9fcf2
Image diff could not be displayed: it is too large. Options to address this: view the blob.
transfert-temps-min.ipynb
+ 39
− 10
View file @ 0ec9fcf2
%% Cell type:markdown id: tags:
# Transfert orbital en temps minimal
* Pour la maison de Fermat.
* Auteur : Olivier Cots
* Date : 05/06/2022
----
On considère le problème de transfert orbital en temps minimal d'un satellite partant d'une orbite initiale elliptique et devant rejoindre une orbite finale circulaire.
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) \\
\ \ \dot{x}_{2}(t) = ~ x_{4}(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]
\ \ 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), \\
\end{array}
\right.
$$
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 :
$$
\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:
## Préliminaires
%% Cell type:markdown id: tags:
### Packages
%% Cell type:code id: tags:
``` julia
using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra
using Plots.PlotMeasures
using SparseArrays
using Printf
int(x) = floor(Int, x);
```
%% Cell type:markdown id: tags:
### Les constantes du problème
%% Cell type:code id: tags:
``` julia
x0 = [-42272.67, 0, 0, -5796.72] # état initial
μ = 5.1658620912*1e12
rf = 42165.0 ;
F_max = 20.0 #100.0
γ_max = F_max*3600.0^2/(2000.0*10^3)
t0 = 0.0
rf3 = rf^3 ;
α = sqrt(μ/rf3);
```
%% Cell type:code id: tags:
``` julia
rayon_Terre = 6371;
sc_sat = 1000;
```
%% Cell type:markdown id: tags:
### Flots de champs de vecteurs
%% Cell type:code id: tags:
``` julia
# Hamiltonian Vector field
mutable struct HV
fun :: Function
end
# Fonction permettant de calculer le flot d'un système hamiltonien
function Flow(hv::HV)
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[:] = hv.fun(z[1:n], z[n+1:2*n])
end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.05)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
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)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end;
# Vector field
mutable struct VF
fun :: Function
end
# Fonction permettant de calculer le flot d'un champ de vecteurs
function Flow(vf::VF)
function rhs!(dx, x, dummy, t)
dx[:] = vf.fun(x)
end
function f(tspan, x0; abstol=1e-12, reltol=1e-12, saveat=0.05)
ode = ODEProblem(rhs!, x0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0, x0, t; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, t), x0; abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol[1:n, end]
end
return f
end;
```
%% Cell type:markdown id: tags:
## Fonctions d'affichage
%% Cell type:markdown id: tags:
### Satellite
%% Cell type:code id: tags:
``` 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
# flamme
@userplot Flamme
@recipe function f(cp::Flamme)
x, y = cp.args
seriestype --> :shape
fillcolor --> :darkorange1
linecolor --> false
x, y
end
function satellite!(pl; subplot=1, thrust=false, 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
# Param bras du satellite
lb = scale.*3*cos(α) # longueur des bras
eb = scale.*cos(α)/30 # demi largeur des bras
xb = 0.0
yb = scale.*sin(α)
# Paramètres corps du satellite
t = range(-α, α, length = 50)
x = scale.*cos.(t);
y = scale.*sin.(t);
Δx = scale.*cos(α)
Δy = scale.*sin(α)
# Paramètres flamme
hF = yb # petite hauteur
HF = 2*hF # grande hauteur
LF = 5.5*Δx # longueur
# Dessin bras du satellite
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,:], subplot=subplot)
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,:], subplot=subplot)
# Dessin corps du satellite
M = T(Rθ*[ x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
M = T(Rθ*[-x'; y'], O); satbody!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
M = T(Rθ*[scale.*[-cos(α), cos(α), cos(α), -cos(α)]'
scale.*[-sin(α), -sin(α), sin(α), sin(α)]'], O);
satbody!(pl, M[1,:], M[2,:], subplot=subplot) # interieur
M = T(Rθ*[ x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord droit
M = T(Rθ*[-x'; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche
M = T(Rθ*[ x'.-2*Δx; y'], O); satline!(pl, M[1,:], M[2,:], subplot=subplot) # bord gauche (droite)
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:], subplot=subplot) # haut
M = T(Rθ*[ scale.*[-cos(α), cos(α)]';
-scale.*[ sin(α), sin(α)]'], O);
satline!(pl, M[1,:], M[2,:], subplot=subplot) # bas
# 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,:], subplot=subplot)
pa2 = T(R(β-π/2)*panneau, v0+v2+vy); pa = T(Rθ*pa2, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa3 = T(R(β-π/2)*panneau, v0+v3+vy); pa = T(Rθ*pa3, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa4 = SA(pa1, β, [xb; yb]); pa = T(Rθ*pa4, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa5 = SA(pa2, β, [xb; yb]); pa = T(Rθ*pa5, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa6 = SA(pa3, β, [xb; yb]); pa = T(Rθ*pa6, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa7 = SA(pa1, 0, [0; 0]); pa = T(Rθ*pa7, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa8 = SA(pa2, 0, [0; 0]); pa = T(Rθ*pa8, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa9 = SA(pa3, 0, [0; 0]); pa = T(Rθ*pa9, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa10 = SA(pa7, -β, [xb; -yb]); pa = T(Rθ*pa10, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa11 = SA(pa8, -β, [xb; -yb]); pa = T(Rθ*pa11, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
pa12 = SA(pa9, -β, [xb; -yb]); pa = T(Rθ*pa12, O); panbody!(pl, pa[1,:], pa[2,:], subplot=subplot)
# Dessin flamme
if thrust
origin = T(Rθ*[Δx, 0.0], O)
thrust = Rθ*[8000, 0.0]
quiver!(pl, [origin[1]], [origin[2]], color=:red,
quiver=([thrust[1]], [thrust[2]]), linewidth=1, subplot=subplot)
end
end;
```
%% Cell type:markdown id: tags:
### Etoiles
%% Cell type:code id: tags:
``` julia
@userplot StarsPlot
@recipe function f(cp::StarsPlot)
xo, yo, Δx, Δy, I, A, m, n = cp.args
seriestype --> :scatter
seriescolor --> :white
seriesalpha --> A
markerstrokewidth --> 0
markersize --> 2*I[3]
xo.+I[2].*Δx./n, yo.+I[1].*Δy./m
end
mutable struct Stars
s
i
m
n
a
end
# Etoiles
function stars(ρ)
n = 200 # nombre de colonnes
m = int(ρ*n) # nombre de lignes
d = 0.03
s = sprand(m, n, d)
i = findnz(s)
s = dropzeros(s)
# alpha
amin = 0.6
amax = 1.0
Δa = amax-amin
a = amin.+rand(length(i[3])).*Δa
return Stars(s, i, m, n, a)
end;
```
%% Cell type:markdown id: tags:
### Trajectoire
%% Cell type:code id: tags:
``` julia
@userplot TrajectoryPlot
@recipe function f(cp::TrajectoryPlot)
t, x, y, tf = cp.args
n = argmin(abs.(t.-tf))
inds = 1:n
seriescolor --> :white
linewidth --> range(0, 3, length = n)
seriesalpha --> range(0, 1.0, length = n)
aspect_ratio --> 1
label --> false
x[inds], y[inds]
end
```
%% Cell type:markdown id: tags:
### Décor
%% Cell type:code id: tags:
``` julia
# Décor
function background(w, h, xmin, xmax, ymin, ymax,
X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, opi, opf)
cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, opi, opf, times, time_current)
# Fond
px = plot(background_color=:gray16, legend = false, aspect_ratio=:equal,
px = plot(background_color=:gray8, legend = false, aspect_ratio=:equal,
size = (w, h), framestyle = :none,
left_margin=-25mm, bottom_margin=-10mm, top_margin=-10mm, right_margin=-10mm,
xlims=(xmin, xmax), ylims=(ymin, ymax) #,
#camera=(0, 90), xlabel = "x", ylabel= "y", zlabel = "z", right_margin=25mm
)
# Etoiles
starsplot!(px, xmin, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin-Δx, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin, ymin-Δy, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin-Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin+Δx, ymin, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin-Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
starsplot!(px, xmin+Δx, ymin+Δy, Δx, Δy, S.i, S.a, S.m, S.n)
# Orbite initiale
nn = length(X2_orb_init)
plot!(px, cx.+X1_orb_init, cy.+X2_orb_init, color = :olivedrab1, linewidth=2, alpha=opi)
# Orbite finale
nn = length(X2_orb_init)
plot!(px, cx.+X1_orb_arr, cy.+X2_orb_arr, color = :turquoise1, linewidth=2, alpha=opf)
# Terre
θ = range(0, 2π, 100)
rT = 6371 #- 1000
rT = rayon_Terre #- 1000
xT = cx .+ rT .* cos.(θ)
yT = cy .+ rT .* sin.(θ)
nn = length(xT)
plot!(px, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
# Soleil
i_current = argmin(abs.(times.-time_current))
e = π/6
β = range(π/4+e, π/4-e, length(times))
ρ = sqrt((0.8*xmax-cx)^2+(0.8*ymax-cy)^2)
cxS = cx + ρ * cos(β[i_current])
cyS = cy + ρ * sin(β[i_current])
θ = range(0, 2π, 100)
rS = 2000
xS = 0.8*xmax .+ rS .* cos.(θ)
yS = 0.8*ymax .+ rS .* sin.(θ)
xS = cxS .+ rS .* cos.(θ)
yS = cyS .+ rS .* sin.(θ)
plot!(px, xS, yS, color = :gold, seriestype=:shape, linewidth=0)
# Point de départ
#plot!(px, [cx+x0[1]], [cy+x0[2]], seriestype=:scatter, color = :white, markerstrokewidth=0)
# Cadre
#plot!(px, [xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], color=:white)
# Angle
#plot!(px, camera=(30,90))
return px
end;
```
%% Cell type:code id: tags:
``` julia
function panneau_control!(px, time_current, times, u, F_max, γ_max, subplot_current)
# We set the (optional) position relative to top-left of the 1st subplot.
# The call is `bbox(x, y, width, height, origin...)`,
# where numbers are treated as "percent of parent".
Δt_control = 2 #0.1*times[end]
xx = 0.08
yy = 0.1
plot!(px, times, F_max.*u[1,:]./γ_max,
inset = (1, bbox(xx, yy, 0.2, 0.15, :top, :left)),
subplot = subplot_current,
bg_inside = nothing,
color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
xguidefontsize=14, legendfontsize=14,
ylims=(-F_max,F_max),
xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
ylabel="u₁ [N]", label=""
)
plot!(px,
[time_current, time_current], [-γ_max,γ_max],
subplot = subplot_current, label="")
plot!(px, times, F_max.*u[2,:]./γ_max,
inset = (1, bbox(xx, yy+0.2, 0.2, 0.15, :top, :left)),
#ticks = nothing,
subplot = subplot_current+1,
bg_inside = nothing,
color=:red, legend=:true, yguidefontrotation=-90, yguidefonthalign = :right,
xguidefontsize=14, legendfontsize=14,
ylims=(-F_max,F_max),
xlims=(time_current-Δt_control, time_current+1.0*Δt_control),
ylabel="u₂ [N]", xlabel="temps [h]", label=""
)
plot!(px,
[time_current, time_current], [-γ_max,γ_max],
subplot = subplot_current+1, label="")
return subplot_current+2
end;
function panneau_information!(px, subplot_current, time_current, i_current,
x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr)
# panneaux information
xx = 0.06
yy = 0.3
plot!(px,
inset = (1, bbox(xx, yy, 0.2, 0.15, :bottom, :left)),
subplot = subplot_current,
bg_inside = nothing, legend=:false,
framestyle=:none
)
a = 0.0
b = 0.0
Δtxt = 0.3
txt = @sprintf("Temps depuis départ [h] = %3.2f", time_current)
annotate!(px, subplot = subplot_current, a, b+2*Δtxt, txt,
annotationcolor=:white, annotationfontsize=14,
annotationhalign=:left)
txt = @sprintf("Vitesse satellite [km/h] = %5.2f", sqrt(v1[i_current]^2+v2[i_current]^2))
annotate!(px, subplot = subplot_current, a, b+Δtxt, txt,
annotationcolor=:white, annotationfontsize=14,
annotationhalign=:left)
txt = @sprintf("Distance à la Terre [km] = %5.2f", sqrt(x1[i_current]^2+x2[i_current]^2))
txt = @sprintf("Distance à la Terre [km] = %5.2f", sqrt(x1[i_current]^2+x2[i_current]^2)-rayon_Terre)
annotate!(px, subplot = subplot_current, a, b, txt,
annotationcolor=:white, annotationfontsize=14,
annotationhalign=:left)
txt = @sprintf("Poussée maximale [N] = %5.2f", F_max)
annotate!(px, subplot = subplot_current, a, b-2*Δtxt, txt,
annotationcolor=:white, annotationfontsize=14,
annotationhalign=:left)
txt = @sprintf("Durée du transfert [h] = %5.2f", tf_transfert)
annotate!(px, subplot = subplot_current, a, b-3*Δtxt, txt,
annotationcolor=:white, annotationfontsize=14,
annotationhalign=:left)
subplot_current = subplot_current+1
plot!(px,
inset = (1, bbox(0.3, 0.03, 0.5, 0.1, :top, :left)),
subplot = subplot_current,
bg_inside = nothing, legend=:false, aspect_ratio=:equal,
xlims=(-4, 4),
framestyle=:none
)
s1 = "Transfert orbital du satellite autour de la Terre\n"
s2 = "de l'orbite elliptique vers l'orbite circulaire\n"
s3 = "en temps minimal"
txt = s1 * s2 * s3
annotate!(px, subplot = subplot_current, 0, 0, txt,
annotationcolor=:white, annotationfontsize=18,
annotationhalign=:center)
rmax = 3*maximum(sqrt.(X1_orb_init.^2+X2_orb_init.^2))
plot!(px, subplot = subplot_current,
0.04.+X1_orb_init./rmax, X2_orb_init./rmax.-0.03,
color = :olivedrab1, linewidth=2)
rmax = 7*maximum(sqrt.(X1_orb_arr.^2+X2_orb_arr.^2))
plot!(px, subplot = subplot_current,
3.3.+X1_orb_arr./rmax, X2_orb_arr./rmax.-0.03,
color = :turquoise1, linewidth=2)
satellite!(px, subplot=subplot_current,
position=[0.67; 0.28], scale=0.08,
rotate=π/2)
# Terre
θ = range(0, 2π, 100)
rT = 0.1
xT = 3.55 .+ rT .* cos.(θ)
yT = 0.28 .+ rT .* sin.(θ)
plot!(px, subplot=subplot_current, xT, yT, color = :dodgerblue1, seriestype=:shape, linewidth=0)
subplot_current = subplot_current+1
xx = 0.0
yy = 0.02
plot!(px,
inset = (1, bbox(xx, yy, 0.12, 0.05, :bottom, :right)),
subplot = subplot_current,
bg_inside = nothing, legend=:false,
framestyle=:none
)
s1 = "Réalisation : Olivier Cots (2022)"
txt = s1
annotate!(px, subplot = subplot_current, 0, 0, txt,
annotationcolor=:gray, annotationfontsize=8, annotationhalign=:left)
return subplot_current+1
end;
```
%% Cell type:markdown id: tags:
### Trajectoire du satellite
%% Cell type:code id: tags:
``` julia
function trajectoire!(px, times, x1, x2, θ, thrust, t_current, cx, cy, pt)
i_current = argmin(abs.(times.-t_current))
# Trajectoire
if t_current>0 && pt[i_current]
trajectoryplot!(px, times, cx.+x1, cy.+x2, t_current)
end
# Satellite
satellite!(px, position=[cx+x1[i_current]; cy+x2[i_current]], scale=sc_sat,
rotate=θ[i_current], thrust=thrust[i_current])
end;
```
%% Cell type:markdown id: tags:
## Animation
%% Cell type:code id: tags:
``` julia
function animation(pre_transfert_data, transfert_data, post_transfert_data; fps=10, nFrame=200)
# pré-transfert
t_pre_transfert = pre_transfert_data.duration
x_pre_transfert = pre_transfert_data.initial_point
# transfert
p0_transfert = transfert_data.initial_adjoint
tf_transfert = transfert_data.duration
# post-trasfert
t_post_transfert = post_transfert_data.duration
# On calcule les orbites initiale et finale
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 = 151
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)
for i in 1:n_theta
theta = Theta[i]
r_orb = p_orb/(1+e_ellipse*cos(theta));
# Orbite initiale
X1_orb_init[i] = r_orb*cos(theta);
X2_orb_init[i] = r_orb*sin(theta);
# Orbite finale
X1_orb_arr[i] = rf*cos(theta) ;
X2_orb_arr[i] = rf*sin(theta);
end
# Centre de la fenêtre
cx = 40000
cy = -7000
# Taille de la fenêtre
w = 1600
h = 900
ρ = h/w
# Limites de la fenêtre
ee = 0.5
xmin = minimum([X1_orb_init; X1_orb_arr]); xmin = xmin - ee * abs(xmin);
xmax = maximum([X1_orb_init; X1_orb_arr]); xmax = xmax + ee * abs(xmax);
ymin = minimum([X2_orb_init; X2_orb_arr]); ymin = ymin - ee * abs(ymin);
ymax = maximum([X2_orb_init; X2_orb_arr]); ymax = ymax + ee * abs(ymax);
Δy = ymax-ymin
Δx = Δy/ρ
Δn = Δx - (xmax-xmin)
xmin = xmin - Δn/2
xmax = xmax + Δn/2
# Trajectoire pré-transfert
traj_pre_transfert = f0((0.0, t_pre_transfert), x_pre_transfert)
times_pre_transfert = traj_pre_transfert.t
n_pre_transfert = size(times_pre_transfert, 1)
x1_pre_transfert = [traj_pre_transfert[1, j] for j in 1:n_pre_transfert ]
x2_pre_transfert = [traj_pre_transfert[2, j] for j in 1:n_pre_transfert ]
v1_pre_transfert = [traj_pre_transfert[3, j] for j in 1:n_pre_transfert ]
v2_pre_transfert = [traj_pre_transfert[4, j] for j in 1:n_pre_transfert ]
# Trajectoire pendant le transfert
traj_transfert = f((t0, tf_transfert), x0, p0_transfert)
times_transfert = traj_transfert.t
n_transfert = size(times_transfert, 1)
x1_transfert = [traj_transfert[1, j] for j in 1:n_transfert ]
x2_transfert = [traj_transfert[2, j] for j in 1:n_transfert ]
v1_transfert = [traj_transfert[3, j] for j in 1:n_transfert ]
v2_transfert = [traj_transfert[4, j] for j in 1:n_transfert ]
u_transfert = zeros(2, length(times_transfert))
for j in 1:n_transfert
u_transfert[:,j] = control(traj_transfert[5:8, j])
end
# post-transfert
x_post_transfert = traj_transfert[:,end]
# Trajectoire post-transfert
traj_post_transfert = f0((0.0, t_post_transfert), x_post_transfert)
times_post_transfert = traj_post_transfert.t
n_post_transfert = size(times_post_transfert, 1)
x1_post_transfert = [traj_post_transfert[1, j] for j in 1:n_post_transfert ]
x2_post_transfert = [traj_post_transfert[2, j] for j in 1:n_post_transfert ]
v1_post_transfert = [traj_post_transfert[3, j] for j in 1:n_post_transfert ]
v2_post_transfert = [traj_post_transfert[4, j] for j in 1:n_post_transfert ]
# Angles de rotation du satellite pendant le pré-transfert
θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])
θ_pre_transfert = range(π/2, mod(θ0, 2*π), n_pre_transfert)
# Angles de rotation du satellite pendant le transfert
θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])
# Angles de rotation du satellite pendant le post-transfert
θ1 = atan(u_transfert[2, end], u_transfert[1, end])
θ2 = atan(-x2_post_transfert[end], -x1_post_transfert[end])
θ_post_transfert = range(mod(θ1, 2*π), mod(θ2, 2*π), n_post_transfert)
# Etoiles
Δx = xmax-xmin
Δy = ymax-ymin
ρ = Δy/Δx
S = stars(ρ)
# nombre total d'images
nFrame = min(nFrame, n_pre_transfert+n_transfert+n_post_transfert);
# Pour l'affichage de la trajectoire globale
times = [times_pre_transfert[1:end-1];
times_pre_transfert[end].+times_transfert[1:end-1];
times_pre_transfert[end].+times_transfert[end].+times_post_transfert[1:end]]
x1 = [x1_pre_transfert[1:end-1]; x1_transfert[1:end-1]; x1_post_transfert[:]]
x2 = [x2_pre_transfert[1:end-1]; x2_transfert[1:end-1]; x2_post_transfert[:]]
v1 = [v1_pre_transfert[1:end-1]; v1_transfert[1:end-1]; v1_post_transfert[:]]
v2 = [v2_pre_transfert[1:end-1]; v2_transfert[1:end-1]; v2_post_transfert[:]]
θ = [ θ_pre_transfert[1:end-1]; θ_transfert[1:end-1]; θ_post_transfert[:]]
# plot thrust on/off
th = [BitArray(zeros(n_pre_transfert-1));
BitArray(ones(n_transfert-1));
BitArray(zeros(n_post_transfert))]
# plot trajectory
pt = [BitArray(ones(n_pre_transfert-1));
BitArray(ones(n_transfert-1));
BitArray(zeros(n_post_transfert))]
# Contrôle sur la trajectoire totale
u_total = hcat([zeros(2, n_pre_transfert-1),
u_transfert[:, 1:n_transfert-1],
zeros(2, n_post_transfert)]...)
# temps total
temps_transfert_global = times[end]
# pas de temps pour le transfert global
if nFrame>1
Δt = temps_transfert_global/(nFrame-1)
else
Δt = 0.0
end
# opacités des orbites initiale et finale
op_initi = [range(0.0, 1.0, length = n_pre_transfert-1);
range(1.0, 0.0, length = n_transfert-1);
zeros(n_post_transfert)]
op_final = [zeros(n_pre_transfert-1);
range(0.0, 1.0, length = n_transfert-1);
range(1.0, 0.0, length = int(n_post_transfert/4));
zeros(n_post_transfert-int(n_post_transfert/4))]
# animation
anim = @animate for i ∈ 1:nFrame
# Δt : pas de temps
# time_current : temps courant de la mission totale à l'itération i
# i_current : indice tel que times[i_current] = time_current
# w, h : width, height de la fenêtre
# xmin, xmax, ymin, ymax : limites des axes du plot principal
# X1_orb_init, X2_orb_init : coordonnées de l'orbite initiale
# X1_orb_arr, X2_orb_arr : coordonnées de l'orbite finale
# cx, cy : coordonées du centre de l'affichage du tranfert
# S : data pour les étoiles
# Δx, Δy : xmax-xmin, ymax-ymin
# times : tous les temps de la mission complète, ie pre-transfert, transfert et post-transfert
# x1, x2 : vecteur de positions du satellite
# θ : vecteur d'orientations du satellite
# th : vecteur de booléens - thrust on/off
# u_total : vecteur de contrôles pour toute la mission
# F_max, γ_max : poussée max
# subplot_current : valeur du subplot courant
# cam_x, cam_y : position de la caméra
# cam_zoom : zoom de la caméra
cam_x = cx
cam_y = cy
cam_zoom = 1
time_current = (i-1)*Δt
i_current = argmin(abs.(times.-time_current))
px = background(w, h, xmin, xmax, ymin, ymax,
X1_orb_init, X2_orb_init, X1_orb_arr, X2_orb_arr,
cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom, op_initi[i_current], op_final[i_current])
cx, cy, S, Δx, Δy, cam_x, cam_y, cam_zoom,
op_initi[i_current], op_final[i_current], times, time_current)
trajectoire!(px, times, x1, x2, θ, th, time_current, cx, cy, pt)
subplot_current = 2
subplot_current = panneau_control!(px, time_current, times, u_total,
F_max, γ_max, subplot_current)
subplot_current = panneau_information!(px, subplot_current, time_current,
i_current, x1, x2, v1, v2, θ, F_max, tf_transfert, X1_orb_init,
X2_orb_init, X1_orb_arr, X2_orb_arr)
end
# enregistrement
gif(anim, "transfert-temps-min.mp4", fps=fps);
nom = "transfert-temps-min-original.mp4"
gif(anim, nom, fps=fps);
gif(anim, "transfert-temps-min.gif", fps=fps);
#return nom
end;
```
%% Cell type:markdown id: tags:
## Calculs
%% Cell type:markdown id: tags:
### Pré-transfert
%% Cell type:code id: tags:
``` julia
# Système non contrôlé : les lois de Kepler
function kepler(x)
n = size(x, 1)
dx = zeros(eltype(x), n)
x1 = x[1]
x2 = x[2]
x3 = x[3]
x4 = x[4]
dsat = norm(x[1:2]); r3 = dsat^3;
dx[1] = x3
dx[2] = x4
dx[3] = -μ*x1/r3
dx[4] = -μ*x2/r3
return dx
end
f0 = Flow(VF(kepler));
```
%% Cell type:code id: tags:
``` julia
# On cherche le point le plus proche de la Terre sur l'orbite initiale
# Ce sera le point de départ de la pré-mission
# On cherche aussi le temps
# données pour la trajectoire pré-transfert
mutable struct PreTransfert
duration
initial_point
end
# Fonction de tir
function depart(xp, tp)
# Integration
x0_ = f0(0.0, xp, tp)
# Conditions
s = zeros(eltype(xp), 5)
s[1] = xp[2]
s[2] = x0_[1] - x0[1]
s[3] = x0_[2] - x0[2]
s[4] = x0_[3] - x0[3]
s[5] = x0_[4] - x0[4]
return s
end;
# Itéré initial
y_guess = [6737.87, -87.10, 272.51, 36364.45, 5.3]
# Jacobienne de la fonction de tir
foo(y) = depart(y[1:4], y[5])
#jfoo(y) = ForwardDiff.jacobian(foo, y)
# Résolution de shoot(p0, tf) = 0
nl_sol = nlsolve(foo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true);
# On récupère la solution si convergence
if converged(nl_sol)
x_pre_transfert = x0 #nl_sol.zero[1:4]
t_pre_transfert = 2.0*nl_sol.zero[5]
pre_transfert_data = PreTransfert(t_pre_transfert, x_pre_transfert)
println("\nFinal solution:\n", nl_sol.zero)
else
error("Not converged")
end
```
%% Cell type:markdown id: tags:
### Transfert
%% Cell type:markdown id: tags:
Le pseudo-hamiltonien est donné par
$$
H(x, p, u) = p_{1}x_3 + p_{2}x_4 + p_{3}\left(-\dfrac{\mu\, x_{1}}{r^{3}} + u_{1}\right) + p_{4}\left(-\dfrac{\mu\, x_{2}}{r^{3}} + u_{2}\right)
$$
et le pseudo système hamiltonien s'écrit
$$
\vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-\frac{\partial H}{\partial x}(x, p, u) \right)
$$
avec
$$
\begin{aligned}
\frac{\partial H}{\partial p}(x, p, u) &= \left( x_3, x_4, -\dfrac{\mu\, x_{1}}{r^{3}} + u_{1}, -\dfrac{\mu\, x_{2}}{r^{3}} + u_{2}\right) \\
\frac{\partial H}{\partial x}(x, p, u) &= \left(\left(\frac{\mu}{r^3} - \frac{3\mu x_1^2}{r^5}\right)p_3 + \frac{3\mu x_1x_2}{r^5}p_4, \frac{3\mu x_1x_2}{r^5}p_3 + \left(\frac{\mu}{r^3} + \frac{3\mu x_2^2}{r^5}\right)p_4, p_1, p_2\right)
\end{aligned}
$$
La maximisation du hamiltonien nous permet d'obtenir le contrôle maximisant
(on suppose que $(p_3,p_4) \neq (0,0)$) :
$$
u(x(t),p(t)) = \frac{\gamma_\mathrm{max}}{\sqrt{p_3^2(t)+p_4^2(t)}}(p_3(t),p_4(t)).
$$
%% Cell type:code id: tags:
``` julia
# Contrôle maximisant
function control(p)
u = zeros(eltype(p),2)
u[1] = p[3]*γ_max/norm(p[3:4])
u[2] = p[4]*γ_max/norm(p[3:4])
return u
end;
# Hamiltonien maximisé
function hfun(x, p)
u = control(p)
nor = norm(x[1:2]); nor3 = nor^3
h = p[1]*x[3] + p[2]*x[4] + p[3]*(-μ*x[1]/nor3 + u[1]) + p[4]*(-μ*x[2]/nor^3 + u[2])
return h
end;
# Système hamiltonien
function hv(x, p)
n = size(x, 1)
hv = zeros(eltype(x), 2*n)
x_1 = x[1]; x_2 = x[2]; x_3 = x[3]; x_4 = x[4]
p_1 = p[1]; p_2 = p[2]; p_3 = p[3]; p_4 = p[4]
norme = norm(x[1:2]); r3 = norme^3; r5 = norme^5
u = control(p)
hv[1] = x_3;
hv[2] = x_4
hv[3] = -μ*x_1/r3 + u[1]
hv[4] = -μ*x_2/r3 + u[2]
hv[5] = μ*((1/r3 - (3*x_1^2)/r5)*p_3 - ((3*x_1*x_2)/r5)*p_4)
hv[6] = μ*(-((3*x_1*x_2)/r5)*p_3 + (1/r3 - (3*x_2^2)/r5)*p_4)
hv[7] = -p_1
hv[8] = -p_2
return hv
end;
f = Flow(HV(hv));
```
%% 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$
donné par
$$
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
$$
H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(cas normal)}
$$
De plus, la condition de transversalité
$$
p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
$$
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.
$$
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{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 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
$z=(x, p)$.
%% Cell type:code id: tags:
``` julia
# Données pour la trajectoire durant le transfert
mutable struct Transfert
duration
initial_adjoint
end
# Fonction de tir
function shoot(p0, tf)
# Integration
xf, pf = f(t0, x0, p0, tf)
# Conditions
s = zeros(eltype(p0), 5)
s[1] = norm(xf[1:2]) - rf
s[2] = xf[3] + α*xf[2]
s[3] = xf[4] - α*xf[1]
s[4] = xf[2]*(pf[1] + α*pf[4]) - xf[1]*(pf[2] - α*pf[3])
s[5] = hfun(xf, pf) - 1
return s
end;
# 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 = [-0.0013615, -7.34989e-6, -5.359923e-5, -0.00858271, 59.8551668] # pour F_max = 20N
# 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 = nl_sol.zero[1:4]
tf_sol = nl_sol.zero[5]
transfert_data = Transfert(tf_sol, p0_sol)
println("\nFinal solution:\n", nl_sol.zero)
else
error("Not converged")
end
```
%% Cell type:markdown id: tags:
### Post-transfert
%% Cell type:code id: tags:
``` julia
# données pour la trajectoire post-transfert
mutable struct PostTransfert
duration
end
t_post_transfert = 20.0
post_transfert_data = PostTransfert(t_post_transfert);
```
%% Cell type:markdown id: tags:
## Création de l'animation
%% Cell type:code id: tags:
``` julia
#nFrame = 100; fps = 10
#nFrame = 10; fps = 1
nFrame = 2000; fps = 20
animation(pre_transfert_data, transfert_data, post_transfert_data, nFrame=nFrame, fps=fps)
```
%% Cell type:code id: tags:
``` julia
;ffmpeg -y -i transfert-temps-min.mp4 -vf format=yuv420p transfert-temps-min-encode.mp4
convert = `ffmpeg -y -i transfert-temps-min-original.mp4 -vf format=yuv420p transfert-temps-min.mp4`
if nFrame > 1
run(convert)
end
```
%% Cell type:markdown id: tags:
[video](./transfert-temps-min.mp4)
%% Cell type:markdown id: tags:
## Annexes
%% Cell type:code id: tags:
``` julia
pl = plot(legend = false, framestyle = :none, aspect_ratio=:equal, size = (400, 400))
scsa = 2000
t = range(0, 2*π, length = 150)
r = 3*scsa
x = r.*cos.(t);
y = r.*sin.(t);
plot!(pl, x, y)
θ = 1*π/4
satellite!(pl; position=0*[3*scsa;0], scale=scsa, rotate=θ, thrust=true)
```
%% Cell type:code id: tags:
``` julia
```
......
transfert-temps-min.mp4
+ 0
− 0
View file @ 0ec9fcf2
No preview for this file type
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment