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

foo

parent 52e3686c
Branches
No related tags found
No related merge requests found
%% Cell type:markdown id:3a61232f-4ef6-47b7-bf22-321b887d49ae tags:
# Orbital transfer in minimum consumption
* Author: Olivier Cots
* Date: 05/06/2022
* For the "Maison de Fermat"
%% Cell type:markdown id:b94045b4-02f3-4b8a-893d-3a3628a3ccff tags:
We consider the orbital transfer in minimum consumption of a satellite around the
Earth, from an initial elliptic orbit to a circular target orbit.
This optimal control problem can be written in the form
$$
\left\lbrace
\begin{array}{l}
\min J(x, u) = \int_{0}^{t_f} \Vert u(t)\Vert \mathrm{d}t \\[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{a. e}, ~~ u(t) = (u_1(t),u_2(t)), \\
\ \ \dot{m}(t) = -\beta\, \Vert u(t)\Vert,
\\[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}, ~~ m(0) = m_0, \\
\ \ 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.
$$
where $r(t)=\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$.
The chosen units are the kilometer for distances et the hour for the times. We set the following parameters:
$$
\mu=5.1658620912 \times 10^{12} \ \mathrm{km}^{3}.\mathrm{h}^{-2}, \quad r_{f} = 42165 \ \mathrm{km}.
$$
The parameter $\gamma_\mathrm{max}$ depends on the maximal thrust $F_\mathrm{max}$ according to
$$
\gamma_\mathrm{max} = \frac{F_\mathrm{max}\times 3600^2}{m~ [g]}
$$
where $m$ is the mass of the satellite that we fix to $m=2000\ \mathrm{kg}$.
%% Cell type:markdown id:f88fa0fe-a75e-47c0-a9bf-e387d44c6f00 tags:
## Packages
%% Cell type:code id:1f8691a7-47ef-41d2-96e3-267af66292fa tags:
``` julia
using MINPACK # fsolve
using LinearAlgebra # norm
```
%% Cell type:code id:6809e33b-a2a6-4508-840b-dd45ed796607 tags:
``` julia
include("flows.jl");
include("plottings.jl");
include("homotopy.jl");
include("commun/flows.jl");
include("commun/plottings.jl");
include("commun/homotopy.jl");
```
%% Cell type:markdown id:a577ddef-dc7b-45ab-93e0-c9bdbd525c0d tags:
## Constants of the problems
%% Cell type:code id:d43f2934-515b-4710-b43c-db4b490b01d1 tags:
``` julia
m0 = 2000.0
x0 = [-42272.67, 0, 0, -5796.72, m0, 0.0]
F_max_20 = 20.0
F_max_30 = 30.0
F_max_40 = 40.0
F_max_60 = 60.0
tf_min_20 = 59.85516688789379 # minimal time for Fmax = 20 N
tf_min_30 = 39.6082272168355
tf_min_40 = 30.02437294598827
tf_min_60 = 19.548304015626
tf_20 = 1.5*tf_min_20
tf_30 = 1.5*tf_min_30
tf_40 = 1.5*tf_min_40
tf_60 = 1.5*tf_min_60
ε_init_1e0 = 1.0
p0_guess_20_a = [0.00668994517486452, 0.0008069227480003886, 0.005884496585344698, -0.04176783251179745, -0.35861897272689147, 0.0]
p0_guess_20_b = [-0.0009737407120528945, 0.0002585826411057628, 0.0018857179140932638, -0.034991324485475844, -0.1693846147531423, 0.0]
p0_guess_30 = [0.0016326413979351802, 0.0004075978609688682, 0.0029724137024974007, -0.019657641630682045, -0.13546101668892516, 0.0]
##p0_guess_40 = [-0.0011276423693230456, 2.706837078747592e-5, 0.00019739651253129466, -0.009847106396009707, -0.034058101939062856, 0.0]
##p0_guess_60 = [-0.0006645768837626829, 4.1034240005117695e-5, 0.00029924282925730477, -0.006676217128598959, -0.02402015651049661, 0.0]
ε_init_1e_1 = 1e-1
p0_guess_20_1e_1 = [-0.00020083445561831767, 8.824275956250577e-5, 0.0006435116910019646, -0.011854368081064033, -0.056208392675317714, 0.0]
ε_init_1e_2 = 1e-2
p0_guess_20_1e_2 = [-9.979024683936127e-5, 7.24625045489202e-5, 0.000528433931997731, -0.010141334052060913, -0.04650571716662274, 0.0]
ε_init_1e_3 = 1e-3
p0_guess_20_1e_3 = [-0.00010295342528640484, 6.955334470904565e-5, 0.0005072188389996865, -0.010062085459039827, -0.045536464232390994, 0.0]
ε_init_5e_4 = 5e-4
p0_guess_20_5e_4 = [0.0]
nx = size(p0_guess_20_a, 1)
Th(F) = F*3600.0^2/(10^3) # u_max
#u_max = F_max*3600.0^2/(10^3)
#γ_max = F_max*3600.0^2/(m0*10^3)
μ = 398600.47 * 3600^2
rf = 42165.0;
rf3 = rf^3
α = sqrt(μ/rf3)
t0 = 0.0
β = 0.0
tol = 1e-12;
```
%% Cell type:markdown id:1677af2f-d61d-4f4a-92fd-078de30b4137 tags:
## The pseudo-Hamiltonian
%% Cell type:markdown id:260b1ace-2c79-4f1a-bb6e-281a300e7fad tags:
The pseudo-Hamiltonian is
$$
H(x, p, u) = -\Vert u \Vert + p_{1}x_3 + p_{2}x_4 + p_{3}\left(-\dfrac{\mu\, x_{1}}{r^{3}} + \frac{u_{1}}{m}\right) + p_{4}\left(-\dfrac{\mu\, x_{2}}{r^{3}} + \frac{u_{2}}{m}\right) - p_5\, \beta\, \Vert u \Vert
$$
%% Cell type:markdown id:4c110234-1a14-4394-9644-fb2c15fc692e tags:
## Logarithmic barrier
To regularize the problem, we introduce
$$
\int_{0}^{t_f} \left( \Vert u(t)\Vert -
\varepsilon \left( \ln( \Vert u(t)\Vert) + \ln(u_\max-\Vert u(t)\Vert) \right) \right)
\mathrm{d}t
$$
and the pseudo-Hamiltonian becomes
$$
H(x, p, u) = - \Vert u \Vert
+ \varepsilon \left( \ln( \Vert u(t)\Vert) + \ln(1-\Vert u(t)\Vert) \right)
+ p_{1}x_3 + p_{2}x_4 + p_{3}\left(-\dfrac{\mu\, x_{1}}{r^{3}} + \frac{u_{1}}{m}\right) + p_{4}\left(-\dfrac{\mu\, x_{2}}{r^{3}} + \frac{u_{2}}{m}\right) - p_5\, \beta\, \Vert u \Vert
$$
%% Cell type:code id:8ea3b2e4-7b75-4a28-bc12-f5733aa9f4e0 tags:
``` julia
# Hamiltonian
# beta = 0 so we do not take it into account
#
function control(x, p, ε, u_max)
φ = [p[3], p[4]]; = norm(φ)
m = x[5]
ψ = -1.0 + (u_max/m)* + p[6] - p[5]*β*u_max
α = (ψ - 2ε + (ψ^2 + 4ε^2))/(2ψ)
u = (α/)*φ
return u
end
#H(x, p, u, ε, u_max) = (-norm(u) + ε*(log(max(norm(u), 1e-3))+log(max(1.0-norm(u), 1e-3)))
H(x, p, u, ε, u_max) = (-norm(u) + ε*(log(norm(u))+log(1.0-norm(u)))
+ p[1]*x[3]
+ p[2]*x[4]
+ p[3]*(-μ*x[1]/norm(x[1:2])^3+(u_max/x[5])*u[1])
+ p[4]*(-μ*x[2]/norm(x[1:2])^3+(u_max/x[5])*u[2])
- p[5]*β*u_max*norm(u)
+ p[6]*norm(u)); # on ajoute la conso dans le hamiltonien: x[6]
H(x, p, ε, u_max) = H(x, p, control(x, p, ε, u_max), ε, u_max)
# Flow
f = Flow(Hamiltonian(H));
```
%% Cell type:code id:88ebc632-6e82-41b5-8fdb-0a5a85adbef1 tags:
``` julia
#=
# Hamiltonian system
# ne fonctionne pas, je ne sais pas pourquoi !
function hv(x, p, u)
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]; m = x[5]
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(x, p, ε)
hv[1] = x_3
hv[2] = x_4
hv[3] = -μ*x_1/r3 + (u_max/m)*u[1]
hv[4] = -μ*x_2/r3 + (u_max/m)*u[2]
hv[5] = -β*u_max*norm(u)
hv[6] = μ*((1/r3 - (3*x_1^2)/r5)*p_3 - ((3*x_1*x_2)/r5)*p_4)
hv[7] = μ*(-((3*x_1*x_2)/r5)*p_3 + (1/r3 - (3*x_2^2)/r5)*p_4)
hv[8] = -p_1
hv[9] = -p_2
hv[10]= (u_max/m^2)*(u[1]*p_3+u[2]*p_4)
return hv
end;
hv(x, p) = hv(x, p, control(x, p, ε))
#f = Flow(HamiltonianVectorField(hv));
=#
```
%% Cell type:markdown id:e14a81c3-2670-4d83-98f7-d1a5df30764b tags:
## Shooting function
%% Cell type:markdown id:d98bcdce-625a-457b-9ae4-b0fe6a797b1a tags:
We introduce
$
\alpha := \sqrt{\frac{\mu}{r_f^3}}.
$
The final condition can be written $c(x(t_f)) = 0$, with $c \colon \mathbb{R}^4 \to \mathbb{R}^3$
given by
$$
c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \alpha\, x_{2}, ~~ x_{4} - \alpha\, x_{1}).
$$
The transversality condition
$$
p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
$$
leads to
$$
\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.
$$
Considering the limit conditions and the transversality condition,
the shooting function is given by
\begin{equation*}
\begin{array}{rlll}
S \colon & \mathbb{R}^4 & \longrightarrow & \mathbb{R}^4 \\
& p_0 & \longmapsto &
S(p_0) := \begin{pmatrix}
c(x(t_f, x_0, p_0)) \\[0.5em]
\Phi(z(t_f, x_0, p_0))
\end{pmatrix}
\end{array}
\end{equation*}
where $z(t_f, x_0, p_0)$ is the solution at time $t_f$ of the feedback pseudo-Hamiltonian system (that is the
pseudo-Hamiltonian system with the maximizing control), starting at time $t_0$ from $(x_0, p_0)$. We recall that
we note $z=(x, p)$.
%% Cell type:code id:ec6f2c03-8752-4246-8a1b-73dfbe097f30 tags:
``` julia
# Données pour la trajectoire durant le transfert
mutable struct Transfert
minimal_time
duration
initial_adjoint
epsilon
Fmax
end
# Fonction de tir
function shoot(p0, tf, ε, u_max)
# Integration
xf, pf = f(t0, x0, p0, tf, ε, u_max, abstol=tol, reltol=tol)
# Conditions
s = zeros(eltype(p0), 6)
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] = pf[5]
s[6] = pf[6]
return s
end;
ε_init = ε_init_1e_3
p0_first_shoot = p0_guess_20_1e_3
tf_first_shoot = tf_20
Fm_first_shoot = F_max_20
tf_min_first_shoot = tf_min_20
# Initial guess
ξ_guess = p0_first_shoot #+ 1e-4.*(-1.0 .+ (2.0 .* rand(Float64, 6)))
# Solve
S(ξ) = shoot(ξ, tf_first_shoot, ε_init, Th(Fm_first_shoot)) # on fixe les valeurs des paramètres
jS(ξ) = ForwardDiff.jacobian(S, ξ)
S!(s, ξ) = ( s[:] = S(ξ); nothing )
jS!(js, ξ) = ( js[:] = jS(ξ); nothing )
println("Initial value of shooting:\n", S(ξ_guess), "\n\n")
nl_sol = fsolve(S!, jS!, ξ_guess, show_trace=true, tol=1e-8); println(nl_sol)
# Retrieves solution
if nl_sol.converged
p0_sol_first_shoot = nl_sol.x[1:nx]
else
error("Not converged")
end
```
%% Cell type:code id:67733fed-63ee-439e-b750-1da3578342b3 tags:
``` julia
# affichage
ode_sol = f((t0, tf_first_shoot), x0, p0_sol_first_shoot, ε_init, Th(Fm_first_shoot), abstol=tol, reltol=tol)
nn = size(ode_sol.t, 1)
uu = zeros(nn, 1)
nx = size(x0, 1)
for j in 1:nn
uu[j] = norm(control(ode_sol[1:nx, j], ode_sol[nx+1:2nx, j], ε_init, Th(Fm_first_shoot)))
end
uu = uu.*Fm_first_shoot
println("consommation = ", ode_sol[6, end]/tf_min_first_shoot)
pp = plot(ode_sol.t, uu, xlabel = "t", ylabel = "||u||", legend=false)
```
%% Cell type:markdown id:74469c3f-47fb-4ebc-a46c-79f2cac890b9 tags:
## Homotopies
%% Cell type:code id:44bbcf2a-6a94-4bab-9a96-174839a8236e tags:
``` julia
Hom(p0, ε) = shoot(p0, tf_first_shoot, ε, Th(Fm_first_shoot))
println("Initial value of homotopy:\n", Hom(p0_sol_first_shoot, ε_init), "\n\n")
if false
p = Path(Hom);
ε₀ = ε_init; ε₁ = 1e-3; path_ε = p(p0_sol_first_shoot, ε₀, ε₁, abstol=1e-12, reltol=1e-12, ftol_proj=1e-6);
p0_guess_ε = path_ε[1:nx, end]
ε_final = ε₁
F_max_final = Fm_first_shoot
tf_final = tf_first_shoot
tf_min_final = tf_min_first_shoot
else
p0_guess_ε = p0_sol_first_shoot
ε_final = ε_init
F_max_final = Fm_first_shoot
tf_final = tf_first_shoot
tf_min_final = tf_min_first_shoot;
end;
```
%% Cell type:code id:8bdc68dc-2108-44f8-beb2-d93b493003f3 tags:
``` julia
# Initial guess
ξ_guess = p0_guess_ε
# Solve
S(ξ) = shoot(ξ, tf_final, ε_final, Th(F_max_final)) # on fixe les valeurs des paramètres
jS(ξ) = ForwardDiff.jacobian(S, ξ)
S!(s, ξ) = ( s[:] = S(ξ); nothing )
jS!(js, ξ) = ( js[:] = jS(ξ); nothing )
println("Initial value of shooting:\n", S(ξ_guess), "\n\n")
nl_sol = fsolve(S!, jS!, ξ_guess, show_trace=true, tol=1e-8); println(nl_sol)
# Retrieves solution
if nl_sol.converged
p0_sol_final = nl_sol.x[1:nx]
transfert_data = Transfert(tf_min_final, tf_final, p0_sol_final, ε_final, F_max_final);
else
error("Not converged")
end
```
%% Cell type:code id:0619b29c-be8c-46b7-b21d-467521a0203b tags:
``` julia
# affichage
ode_sol = f((t0, tf_final), x0, p0_sol_final, ε_final, Th(F_max_final), abstol=tol, reltol=tol)
nn = size(ode_sol.t, 1)
uu = zeros(nn, 1)
nx = size(x0, 1)
for j in 1:nn
uu[j] = norm(control(ode_sol[1:nx, j], ode_sol[nx+1:2nx, j], ε_final, Th(F_max_final)))
end
uu = uu.*F_max_final
println("consommation = ", ode_sol[6, end]/tf_min_final)
pp = plot(ode_sol.t, uu, xlabel = "t", ylabel = "||u||", legend=false)
```
%% Cell type:code id:d0f3a6df-b0f1-476c-892d-2645e6ebeb7a tags:
``` julia
#=
# Homotopy on tf and Fmax (u_max for shooting)
F_max_init = F_max_60
F_max_final = F_max_40
tf_init = tf_60
tf_final = tf_40
ω(λ) = (1.0-λ).*[tf_init, Th(F_max_init)] + λ.*[tf_final, Th(F_max_final)]
Hom(p0, λ) = shoot(p0, ω(λ)[1], ε_final, ω(λ)[2])
println("Initial value of homotopy:\n", Hom(p0_sol_final, 0.0), "\n\n")
p = Path(Hom);
λ₀ = 0.; λ₁ = 1.; path_Fmax = p(p0_sol_final, λ₀, λ₁);
=#
```
%% Cell type:code id:164db52a-546f-4aef-8d1d-eac4e575be32 tags:
``` julia
#=
p0_guess_F = path_Fmax[1:nx, end]
#ε_final = ε₁
#F_max_final = F_max_40
#tf_final = tf_60
tf_min_final = tf_min_40
# Initial guess
ξ_guess = p0_guess_F
# Solve
S(ξ) = shoot(ξ, tf_final, ε_final, Th(F_max_final)) # on fixe les valeurs des paramètres
jS(ξ) = ForwardDiff.jacobian(S, ξ)
S!(s, ξ) = ( s[:] = S(ξ); nothing )
jS!(js, ξ) = ( js[:] = jS(ξ); nothing )
println("Initial value of shooting:\n", S(ξ_guess), "\n\n")
nl_sol = fsolve(S!, jS!, ξ_guess, show_trace=true, tol=1e-8); println(nl_sol)
# Retrieves solution
if nl_sol.converged
p0_sol_final = nl_sol.x[1:nx]
transfert_data = Transfert(tf_min_final, tf_final, p0_sol_final, ε_final, F_max_final);
else
error("Not converged")
end
=#
```
%% Cell type:code id:37a0beea-c5a6-48e7-a825-a828f1d98300 tags:
``` julia
#=
# affichage
ode_sol = f((t0, tf_final), x0, p0_sol_final, ε_final, Th(F_max_final))
nn = size(ode_sol.t, 1)
uu = zeros(nn, 1)
nx = size(x0, 1)
for j in 1:nn
uu[j] = norm(control(ode_sol[1:nx, j], ode_sol[nx+1:2nx, j], ε_final, Th(F_max_final)))
end
uu = uu.*F_max_final
pp = plot(ode_sol.t, uu, xlabel = "t", ylabel = "||u||", legend=false)
=#
```
%% Cell type:markdown id:3038276d-9a9e-4f3f-b193-06b9b0f13fdd tags:
## Pre and post-transfer
%% Cell type:code id:f62cca9f-2099-409b-a873-13f6041b01f1 tags:
``` julia
# Kepler's law
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(VectorField(kepler));
```
%% Cell type:code id:81f5787b-0609-4b25-8ff1-6fdd8d164a83 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
ξ_guess = [6738.200652584018, 1.3618799608695503e-21, 2.8839710428042565e-7, 36366.2117341436, 5.302395297743802]
# Solve
foo(ξ) = depart(ξ[1:4], ξ[5])
jfoo(ξ) = ForwardDiff.jacobian(foo, ξ)
foo!(s, ξ) = ( s[:] = foo(ξ); nothing )
jfoo!(js, ξ) = ( js[:] = jfoo(ξ); nothing )
nl_sol = fsolve(foo!, jfoo!, ξ_guess, show_trace=true); println(nl_sol)
# Retrieves solution
if nl_sol.converged
x_pre_transfert = x0
t_pre_transfert = 2.0*nl_sol.x[5]
pre_transfert_data = PreTransfert(t_pre_transfert, x_pre_transfert);
else
error("Not converged")
end
```
%% Cell type:code id:d909fca0-a90b-4b63-8f4a-1d32f3b9ab68 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:98f86df9-1add-4b94-9273-741f083d48f2 tags:
## Animation
%% Cell type:code id:c56dc45b-d0e0-427f-944b-239cb0223c31 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
tf_min = transfert_data.minimal_time
ε = transfert_data.epsilon
F_max = transfert_data.Fmax
# 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.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 ]
c_pre_transfert = zeros(1, n_pre_transfert)
# Trajectoire pendant le transfert
traj_transfert = f((t0, tf_transfert), x0, p0_transfert, ε, Th(F_max), abstol=tol, reltol=tol)
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, n_transfert)
n = size(x0, 1)
for j in 1:n_transfert
u_transfert[:,j] = control(traj_transfert[1:n, j], traj_transfert[n+1:2n, j], ε, Th(F_max))
end
#u_transfert = u_transfert # ||u|| ≤ 1
c_transfert = [traj_transfert[6, j] for j in 1:n_transfert ]
# post-transfert
x_post_transfert = traj_transfert[1:4,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 ]
c_post_transfert = zeros(1, n_post_transfert) .+ c_transfert[end]
# Angles de rotation du satellite pendant le pré-transfert
# Et poussée normalisée entre 0 et 1
θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])
θ_pre_transfert = range(π/2, mod(θ0, 2*π), length=n_pre_transfert)
F_pre_transfert = zeros(1, n_pre_transfert)
# Angles de rotation du satellite pendant le transfert
θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])
F_transfert = zeros(1, n_transfert)
th_transfert = BitArray(ones(n_transfert))
for j in 1:n_transfert
F_transfert[j] = norm(u_transfert[:,j]) # in [0, 1]
if F_transfert[j] < 1e-1 # on n'affiche pas la poussée (sur le saltelitte)
# si elle est faible car on est sur un pb régularisé mais on veut afficher du bang-bang
th_transfert[j] = false
end
end
# 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*π), length=n_post_transfert)
F_post_transfert = zeros(1, 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[:]]
F = [ F_pre_transfert[1:end-1]; F_transfert[1:end-1]; F_post_transfert[:]]
cc = [c_pre_transfert[1:end-1]; c_transfert[1:end-1]; c_post_transfert[:]]
# plot thrust on/off
th = [BitArray(zeros(n_pre_transfert-1));
th_transfert[1: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))]
println("Fps = ", fps)
println("NFrame = ", nFrame)
println("Temps transfert = ", temps_transfert_global)
println("Vitesse = ", nFrame/(temps_transfert_global*fps))
# 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], times, time_current)
trajectoire!(px, times, x1, x2, θ, F, th, time_current, cx, cy, pt)
subplot_current = 2
subplot_current = panneau_control!(px, time_current, times, u_total,
F_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, pb=consomin, tf_min=tf_min, conso=cc)
end
# enregistrement
gif(anim, "transfert-conso-min-original.mp4", fps=fps);
gif(anim, "transfert-conso-min.gif", fps=fps);
end;
```
%% Cell type:code id:f0fcd66e-d2f1-46fe-b589-972180b384be tags:
``` julia
function animation(pre_transfert_data, transfert_data, post_transfert_data, vitesse, fps)
t_pre_transfert = pre_transfert_data.duration
t_transfert = transfert_data.duration
t_post_transfert = post_transfert_data.duration
t = t_pre_transfert + t_transfert + t_post_transfert
d = vitesse * t
N = floor(Int, d*fps)
println("Fps = ", fps)
println("NFrame = ", N)
println("Temps transfert = ", t)
println("Vitesse = ", vitesse)
println("")
animation(pre_transfert_data, transfert_data, post_transfert_data, nFrame=N, fps=fps)
end;
```
%% Cell type:code id:40ff10ea-988f-4e81-8f51-a2f1b8ff28cc tags:
``` julia
if false
#nFrame = 1; fps = 1
nFrame = 100; d = 30; fps = floor(Int, nFrame/d)
#nFrame = 2000; d = 90; fps = floor(Int, nFrame/d)
animation(pre_transfert_data, transfert_data, post_transfert_data, nFrame=nFrame, fps=fps)
if nFrame == 2000
convert = `ffmpeg -y -i transfert-conso-min-original.mp4 -vf format=yuv420p transfert-conso-min.mp4`
run(convert)
end
else
vitesse = 0.8421606675612892
fps = 22
animation(pre_transfert_data, transfert_data, post_transfert_data, vitesse, fps)
convert = `ffmpeg -y -i transfert-conso-min-original.mp4 -vf format=yuv420p transfert-conso-min.mp4`
run(convert)
end
```
%% Cell type:code id:7eebd6f1-aedf-43f5-9fd1-ecacce3a5e5a tags:
``` julia
```
......
%% Cell type:markdown id:baa07c4f-8a11-49d2-b052-e64d34594294 tags:
# Orbital transfer in minimum time
* Author: Olivier Cots
* Date: 05/06/2022
* For the "Maison de Fermat"
%% Cell type:markdown id:b32c3c7c-dcc2-4dd9-9143-b4daf05f1c88 tags:
We consider the orbital transfer in minimum time of a satellite around the Earth,
from an initial elliptic orbit to a circular target orbit.
This optimal control problem can be written in the form
$$
\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{a. e}, ~~ 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.
$$
where $r(t)=\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$.
The chosen units are the kilometer for distances et the hour for the times. We set the following parameters:
$$
\mu=5.1658620912 \times 10^{12} \ \mathrm{km}^{3}.\mathrm{h}^{-2}, \quad r_{f} = 42165 \ \mathrm{km}.
$$
The parameter $\gamma_\mathrm{max}$ depends on the maximal thrust $F_\mathrm{max}$ according to
$$
\gamma_\mathrm{max} = \frac{F_\mathrm{max}\times 3600^2}{m ~ [g]}
$$
where $m$ is the mass of the satellite that we fix to $m=2000\ \mathrm{kg}$.
%% Cell type:markdown id:a552de56-944a-412b-97b4-a71c52123786 tags:
## Packages
%% Cell type:code id:104fc457-beaa-4783-8ce2-bb32ff0ab9e3 tags:
``` julia
using MINPACK # fsolve
using LinearAlgebra # norm
```
%% Cell type:code id:0995cf52-227b-4c1f-88cd-ec3f6bcde259 tags:
``` julia
include("../commun/flows.jl");
include("../commun/plottings.jl");
include("commun/flows.jl");
include("commun/plottings.jl");
```
%% Cell type:markdown id:9b5d4e41-f03e-4e62-9970-8c0c24642119 tags:
## Constants of the problems
%% Cell type:code id:f39bc943-a1d2-47a7-81ea-3d766c9d11a6 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:markdown id:80679ed8-77de-43aa-b542-6a5f7733e86b tags:
## The Hamiltonian system
%% Cell type:markdown id:2139105f-a494-4692-8225-46c165e1fae0 tags:
The pseudo-Hamiltonian is
$$
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)
$$
and the pseudo-Hamiltonian system is given by
$$
\vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-\frac{\partial H}{\partial x}(x, p, u) \right)
$$
where
$$
\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}
$$
The maximization condition from the Pontryagin Maximum Principle gives the maximizing control
(we assume that $(p_3,p_4)$ never vanishes) :
$$
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:650273cc-a4fc-46f9-ba02-a65f1a57ed90 tags:
``` julia
# Maximizing control
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;
# Maximized Hamiltonian
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;
# Flow
z = Flow(Hamiltonian(hfun));
# Hamiltonian system
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(HamiltonianVectorField(hv));
```
%% Cell type:markdown id:6bda152c-5c5f-4de3-aad9-6e7b0f7975f4 tags:
## The shooting function
%% Cell type:markdown id:7a35da8c-771e-4e57-8413-fc615d140474 tags:
We introduce
$
\alpha := \sqrt{\frac{\mu}{r_f^3}}.
$
The final condition can be written $c(x(t_f)) = 0$, with $c \colon \mathbb{R}^4 \to \mathbb{R}^3$
given by
$$
c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \alpha\, x_{2}, ~~ x_{4} - \alpha\, x_{1}).
$$
The final time being free, we get the condition at the final time
$$
H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(cas normal)}
$$
Besides, the transversality condition
$$
p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
$$
leads to
$$
\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.
$$
Considering the limit conditions, the final condition on the pseudo-Hamiltonian and the transversality condition,
the shooting function is given by
\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*}
where $z(t_f, x_0, p_0)$ is the solution at time $t_f$ of the feedback pseudo-Hamiltonian system (that is the
pseudo-Hamiltonian system with the maximizing control), starting at time $t_0$ from $(x_0, p_0)$. We recall that
we note $z=(x, p)$.
%% Cell type:code id:4deeb191-35a8-4f5b-9dca-efa1e7b1ac8b 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 = z(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;
# Initial guess
#y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4] # for F_max = 100N
ξ_guess = [-0.0013615, -7.34989e-6, -5.359923e-5, -0.00858271, 59.8551668] # for F_max = 20N
# Solve
foo(ξ) = shoot(ξ[1:4], ξ[5])
jfoo(ξ) = ForwardDiff.jacobian(foo, ξ)
foo!(s, ξ) = ( s[:] = foo(ξ); nothing )
jfoo!(js, ξ) = ( js[:] = jfoo(ξ); nothing )
nl_sol = fsolve(foo!, jfoo!, ξ_guess, show_trace=true); println(nl_sol)
# Retrieves solution
if nl_sol.converged
p0_sol = nl_sol.x[1:4]
tf_sol = nl_sol.x[5]
transfert_data = Transfert(tf_sol, p0_sol);
else
error("Not converged")
end
```
%% Cell type:markdown id:ae9cd7d1-97b2-465f-8f05-f2eeea9b1def tags:
## Pre and post-transfer
%% Cell type:code id:73dcf5ff-4ff2-4894-9407-720bed9b4983 tags:
``` julia
# Kepler's law
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(VectorField(kepler));
```
%% Cell type:code id:04f4e32f-ce8a-4e1a-90a1-b3dc0af2a956 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
ξ_guess = [6738.200652584018, 1.3618799608695503e-21, 2.8839710428042565e-7, 36366.2117341436, 5.302395297743802]
# Solve
foo(ξ) = depart(ξ[1:4], ξ[5])
jfoo(ξ) = ForwardDiff.jacobian(foo, ξ)
foo!(s, ξ) = ( s[:] = foo(ξ); nothing )
jfoo!(js, ξ) = ( js[:] = jfoo(ξ); nothing )
nl_sol = fsolve(foo!, jfoo!, ξ_guess, show_trace=true); println(nl_sol)
# Retrieves solution
if nl_sol.converged
x_pre_transfert = x0
t_pre_transfert = 2.0*nl_sol.x[5]
pre_transfert_data = PreTransfert(t_pre_transfert, x_pre_transfert);
else
error("Not converged")
end
```
%% Cell type:code id:00a00a7d-a161-4327-9ea3-f8e4a34bfb77 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:0686c07f-aea9-405a-95bd-af42f09370b6 tags:
## Animation
%% Cell type:code id:e61c1bee-88b4-4b57-8a09-1249f31978e0 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.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
u_transfert = u_transfert./γ_max # ||u|| ≤ 1
# 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
# Et poussée normalisée entre 0 et 1
θ0 = atan(u_transfert[2, 1], u_transfert[1, 1])
θ_pre_transfert = range(π/2, mod(θ0, 2*π), length=n_pre_transfert)
F_pre_transfert = zeros(1, n_pre_transfert)
# Angles de rotation du satellite pendant le transfert
θ_transfert = atan.(u_transfert[2, :], u_transfert[1, :])
F_transfert = zeros(1, n_transfert)
for j in 1:n_transfert
F_transfert[j] = norm(u_transfert[:,j])
end
# 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*π), length=n_post_transfert)
F_post_transfert = zeros(1, 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[:]]
F = [ F_pre_transfert[1:end-1]; F_transfert[1:end-1]; F_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))]
println("Fps = ", fps)
println("NFrame = ", nFrame)
println("Temps transfert = ", temps_transfert_global)
println("Vitesse = ", nFrame/(temps_transfert_global*fps))
# 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], times, time_current)
trajectoire!(px, times, x1, x2, θ, F, th, time_current, cx, cy, pt)
subplot_current = 2
subplot_current = panneau_control!(px, time_current, times, u_total,
F_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-original.mp4", fps=fps);
gif(anim, "transfert-temps-min.gif", fps=fps);
end;
```
%% Cell type:code id:a23d5c91-72e5-4a93-b689-8227887f0551 tags:
``` julia
#nFrame = 1; fps = 1
#nFrame = 100; d = 30; fps = floor(Int, nFrame/d)
nFrame = 2000; d = 90; fps = floor(Int, nFrame/d)
animation(pre_transfert_data, transfert_data, post_transfert_data, nFrame=nFrame, fps=fps)
```
%% Cell type:code id:acc4155b-8c6a-4e06-a823-8c22fd9ff72a tags:
``` julia
if nFrame == 2000
convert = `ffmpeg -y -i transfert-temps-min-original.mp4 -vf format=yuv420p transfert-temps-min.mp4`
run(convert)
end
```
%% Cell type:code id:c663a64a-1e53-43ba-8fbf-86fb073cc526 tags:
``` julia
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment