OptimalControl.jl
OptimalControl.jl copied to clipboard
control-toolbox
[General] Can't add more than 3 arcs in the direct optimisation
Voici le code.
Le debut est pour intégrer les toolbox nécessaire + définir les paramètres. Ensuite je formule 3 pbs: 1 arc fonctionne, 3 arcs aussi, 4 non...
Tout à la fin il y a une fonction de plot que je définis, pas nécessaire.
using SPICE
using OptimalControl
using NLPModelsIpopt
# using Interpolations
using JLD2
# using OrdinaryDiffEq
using DataInterpolations
using CSV
using JLD2
using Plots
using Downloads: download
# SPICE files that you have to download but only once
const LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
const SPK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp"
const PCK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc"
# Download kernels
download(LSK, "naif0012.tls")
download(SPK, "de440.bsp")
download(PCK, "pck00010.tpc")
furnsh("naif0012.tls") # Leap seconds kernel
furnsh("de430.bsp") # Ephemeris kernel
furnsh("pck00010.tpc") # PCK kernel for planetary constants
# Astrodynamical parameters
G = 6.6743e-11 / 1000^3 # km^3 kg^-1 s^-2
M_moon = 7.346e22 # kg
M_earth = 5.972168e24 # kg
M_sun = 1.989e30 # kg
AU = 1.495978707e11 / 1e3 # km
LU = 384399 # km
TU = sqrt(LU^3 / G / (M_moon + M_earth))
MUnit = M_moon + M_earth
mu_Earth = 3.9860044188e14 / 1000^3 / LU^3 * TU^2
mu_Moon = 4.90486959e12 / 1000^3 / LU^3 * TU^2
mu_Sun = 1.327124400189e20 / 1000^3 / LU^3 * TU^2
mass_moon = M_moon / MUnit
mass_earth = M_earth / MUnit
mass_sun = M_sun / MUnit
# Spacecraft parameters
g0 = 9.80665 / 1000 / LU * TU^2
Isp = 1660 / TU
wet_mass = 344 / MUnit
fuel_mass = 90 / MUnit
dry_mass = (wet_mass - fuel_mass)
wet_by_dry = (dry_mass + fuel_mass) / dry_mass
Tmax = 0.090 / 1000 / MUnit / LU * TU^2 #thrust force in N
Q = Tmax / Isp / g0
mass_scaling_factor = wet_mass
# Mission data that I just copy paste here for simplicity
t0_J200 = 1.107747244e9
tf_J200 = 1.112803048e9
N = 15000
t_values = LinRange(t0_J200, tf_J200, N-3)
t_values = vcat(-5e9, 0.0, collect(t_values), tf_J200 * 5)
states_moon = zeros(6, N)
states_earth = zeros(6, N)
states_sun = zeros(6, N)
t0 = t0_J200 / TU
tf = tf_J200 / TU
# Interpolations of the bodies positions/velocities
for i = 1:N
states_moon[:,i], _ = spkezr("MOON", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
states_earth[:,i], _ = spkezr("EARTH", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
states_sun[:,i], _ = spkezr("SUN", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
end
# Interpolation of the position and velocity of the celestial bodies (Sun, Moon, Earth)
x_moon = DataInterpolations.LinearInterpolation(states_moon[1,:] / LU, t_values / TU)
y_moon = DataInterpolations.LinearInterpolation(states_moon[2,:] / LU, t_values / TU)
z_moon = DataInterpolations.LinearInterpolation(states_moon[3,:] / LU, t_values / TU)
x_earth = DataInterpolations.LinearInterpolation(states_earth[1,:] / LU, t_values / TU)
y_earth = DataInterpolations.LinearInterpolation(states_earth[2,:] / LU, t_values / TU)
z_earth = DataInterpolations.LinearInterpolation(states_earth[3,:] / LU, t_values / TU)
x_sun = DataInterpolations.LinearInterpolation(states_sun[1,:] / LU, t_values / TU)
y_sun = DataInterpolations.LinearInterpolation(states_sun[2,:] / LU, t_values / TU)
z_sun = DataInterpolations.LinearInterpolation(states_sun[3,:] / LU, t_values / TU)
# Times
t0 = 2952.5065962806448
t1 = 2953.955962786101
t2 = 2954.4423792964794
t3 = 2955.398985817685
t4 = 2957.4519685112937
t5 = 2959.0986010622196
t6 = 2960.6895789585565
# States at the previous times
ap0 = [-0.992463207586249; 0.7937735622182638; 0.07902677416947182; 0.4497665242298916; 0.5768992856452368; -0.08990077540870026]
occult1 = [0.024011456548615025; 1.0839463211562548; -0.06721110889958525; 0.8370532723530986; -0.2849379073937811; -0.09055054145691602]
lunar2 = [0.42191676845234877; 0.8365166413216278; -0.10058044976319964; 1.0641134073373073; -0.7420861511380829; 0.2221102307286784]
ap3 = [0.8866648294530993; 0.07703956973358668; 0.28758014926093095; 0.010909379448810359; -1.0104415437867642; 0.3065975859809964]
lunar4 = [-0.5240519542473246; -0.8727170323241018; 0.06983526263136841; -1.0656424593386877; 0.4013541226489389; -0.24177176981236737]
ap5 = [-1.2845625006191574; 0.07452729181167739; -0.02061966816776342; 0.03379116945381267; 0.7236740339189934; -0.05964259399261467]
occult6 = [-0.5063291057829847; 0.9375405520827109; -0.0856392904912281; 0.906825723975573; 0.17758145410610057; -0.009952663265456345]
# Changement de variable pour le temps
function phi(s, t0, t1)
return t0 + (t1 - t0) * s
end
# Dynamique
function Fu(t, x, u)
# Earth & Moon & Sun gravity
r = x[1:3]
v = x[4:6]
m = x[7] * mass_scaling_factor
x_e = [x_earth(t), y_earth(t), z_earth(t)]
x_m = [x_moon(t), y_moon(t), z_moon(t)]
x_s = [x_sun(t), y_sun(t), z_sun(t)]
normrE = sqrt(sum((x_e - r).^2))
normrM = sqrt(sum((x_m - r).^2))
normrS = sqrt(sum((x_s - r).^2))
normrSE = sqrt(sum((x_e - x_s).^2))
normrMS = sqrt(sum((x_m - x_s).^2))
dv = mu_Earth / normrE^3 * (x_e - r) + mu_Moon / normrM^3 * (x_m - r) + mu_Sun / normrS^3 * (x_s - r) - (1 / (mass_earth + mass_moon)) * ( x_s * (mu_Sun * mass_earth / normrSE^3 + mu_Moon * mass_sun / normrMS^3) - x_m * (mu_Moon * mass_sun / normrMS^3) - x_e * (mu_Sun * mass_earth / normrSE^3))
dvu = Tmax * u / m
dm = - Q * sqrt(sum(u.^2)) / mass_scaling_factor
dx = [v; dv + dvu; dm]
return dx
end
# Pb avec 1 arc
@def ocp1arc begin
s ∈ [ 0, 1 ], time
y ∈ R^7, state
w ∈ R^3, control
x1 = y[1:7]
u1 = w[1:3]
mf = y[end]
-5 ≤ x1[1](s) ≤ 5
-5 ≤ x1[2](s) ≤ 5
-5 ≤ x1[3](s) ≤ 5
-5 ≤ x1[4](s) ≤ 5
-5 ≤ x1[5](s) ≤ 5
-5 ≤ x1[6](s) ≤ 5
0.5 ≤ x1[7](s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1
x1[1:6](0) == occult1
x1[7](0) == 1
x1[1:6](1) == lunar2
ẏ(s) == (t2 - t1) * Fu(phi(s, t1, t2), x1(s), u1(s))
mf(1) → max
end
sol = solve(ocp1arc)
plot(sol)
# sol = solve(ocp1arc, grid_size = 100, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5)
# 3 Arcs
@def ocp3arc begin
s ∈ [ 0, 1 ], time
y ∈ R^21, state
w ∈ R^9, control
x1 = y[1:7]
x2 = y[8:14]
x3 = y[15:21]
u1 = w[1:3]
u2 = w[4:6]
u3 = w[7:9]
mf = y[end]
-5 ≤ x1[1](s) ≤ 5
-5 ≤ x1[2](s) ≤ 5
-5 ≤ x1[3](s) ≤ 5
-5 ≤ x1[4](s) ≤ 5
-5 ≤ x1[5](s) ≤ 5
-5 ≤ x1[6](s) ≤ 5
0.5 ≤ x1[7](s) ≤ 1 # careful here ! dry mass > 0
-5 ≤ x2[1](s) ≤ 5
-5 ≤ x2[2](s) ≤ 5
-5 ≤ x2[3](s) ≤ 5
-5 ≤ x2[4](s) ≤ 5
-5 ≤ x2[5](s) ≤ 5
-5 ≤ x2[6](s) ≤ 5
0.5 ≤ x2[7](s) ≤ 1 # careful here ! dry mass > 0
-5 ≤ x3[1](s) ≤ 5
-5 ≤ x3[2](s) ≤ 5
-5 ≤ x3[3](s) ≤ 5
-5 ≤ x3[4](s) ≤ 5
-5 ≤ x3[5](s) ≤ 5
-5 ≤ x3[6](s) ≤ 5
0.5 ≤ x3[7](s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1
0 ≤ sum(u2(s).^2) ≤ 1
0 ≤ sum(u3(s).^2) ≤ 1
x1[1:6](0) == ap0
x1[7](0) == 1
x1[1:6](1) == occult1
x2(0) - x1(1) == [0, 0, 0, 0, 0, 0, 0]
x2[1:6](1) == lunar2
x3(0) - x2(1) == [0, 0, 0, 0, 0, 0, 0]
x3[1:6](1) == ap3
ẏ(s) == [(t1 - t0) * Fu(phi(s, t0, t1), x1(s), u1(s)); (t2 - t1) * Fu(phi(s, t1, t2), x2(s), u2(s)); (t3 - t2) * Fu(phi(s, t2, t3), x3(s), u3(s))]
mf(1) → max
end
sol = solve(ocp3arc, grid_size = 100, tol = 1e-6, disc_method=:gauss_legendre_2)
plot(sol)
# plot_trajNV_multiarc(sol, 3, [t0,t1,t2,t3])
# 4 arcs
@def ocp4arc begin
s ∈ [ 0, 1 ], time
y ∈ R^28, state
w ∈ R^12, control
x1 = y[1:7]
x2 = y[8:14]
x3 = y[15:21]
x4 = y[22:28]
u1 = w[1:3]
u2 = w[4:6]
u3 = w[7:9]
u4 = w[10:12]
mf = y[end]
-5 ≤ x1[1](s) ≤ 5
-5 ≤ x1[2](s) ≤ 5
-5 ≤ x1[3](s) ≤ 5
-5 ≤ x1[4](s) ≤ 5
-5 ≤ x1[5](s) ≤ 5
-5 ≤ x1[6](s) ≤ 5
0.5 ≤ x1[7](s) ≤ 1 # careful here ! dry mass > 0
-5 ≤ x2[1](s) ≤ 5
-5 ≤ x2[2](s) ≤ 5
-5 ≤ x2[3](s) ≤ 5
-5 ≤ x2[4](s) ≤ 5
-5 ≤ x2[5](s) ≤ 5
-5 ≤ x2[6](s) ≤ 5
0.5 ≤ x2[7](s) ≤ 1 # careful here ! dry mass > 0
-5 ≤ x3[1](s) ≤ 5
-5 ≤ x3[2](s) ≤ 5
-5 ≤ x3[3](s) ≤ 5
-5 ≤ x3[4](s) ≤ 5
-5 ≤ x3[5](s) ≤ 5
-5 ≤ x3[6](s) ≤ 5
0.5 ≤ x3[7](s) ≤ 1 # careful here ! dry mass > 0
-5 ≤ x4[1](s) ≤ 5
-5 ≤ x4[2](s) ≤ 5
-5 ≤ x4[3](s) ≤ 5
-5 ≤ x4[4](s) ≤ 5
-5 ≤ x4[5](s) ≤ 5
-5 ≤ x4[6](s) ≤ 5
0.5 ≤ x4[7](s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1
0 ≤ sum(u2(s).^2) ≤ 1
0 ≤ sum(u3(s).^2) ≤ 1
0 ≤ sum(u4(s).^2) ≤ 1
x1[1:6](0) == ap0
x1[7](0) == 1
x1[1:6](1) == occult1
x2(0) - x1(1) == [0, 0, 0, 0, 0, 0, 0]
x2[1:6](1) == lunar2
x3(0) - x2(1) == [0, 0, 0, 0, 0, 0, 0]
x3[1:6](1) == ap3
x4(0) - x3(1) == [0, 0, 0, 0, 0, 0, 0]
x4[1:6](1) == lunar4
ẏ(s) == [(t1 - t0) * Fu(phi(s, t0, t1), x1(s), u1(s)); (t2 - t1) * Fu(phi(s, t1, t2), x2(s), u2(s)); (t3 - t2) * Fu(phi(s, t2, t3), x3(s), u3(s)); (t4 - t3) * Fu(phi(s, t3, t4), x4(s), u4(s))]
mf(1) → max
end
sol = solve(ocp4arc)
function plot_trajNV_multiarc(sol, nb_arcs, tvec)
sol_time = sol.time_grid.value
Nsol = length(sol_time)
sol_state = zeros(7 * nb_arcs, Nsol)
sol_control = zeros(4 * nb_arcs, Nsol)
for i in range(1, length(sol.time_grid.value))
sol_state[:, i] = sol.state.value(sol.time_grid.value[I])
end
x_vals_opti = zeros(nb_arcs, Nsol)
y_vals_opti = zeros(nb_arcs, Nsol)
z_vals_opti = zeros(nb_arcs, Nsol)
vx_vals_opti = zeros(nb_arcs, Nsol)
vy_vals_opti = zeros(nb_arcs, Nsol)
vz_vals_opti = zeros(nb_arcs, Nsol)
m_vals_opti = zeros(nb_arcs, Nsol)
for i in range(1, nb_arcs)
x_vals_opti[i, :] = sol_state[1 + 7 * (i - 1), :]
y_vals_opti[i, :] = sol_state[2 + 7 * (i - 1), :]
z_vals_opti[i, :] = sol_state[3 + 7 * (i - 1), :]
vx_vals_opti[i, :] = sol_state[4 + 7 * (i - 1), :]
vy_vals_opti[i, :] = sol_state[5 + 7 * (i - 1), :]
vz_vals_opti[i, :] = sol_state[6 + 7 * (i - 1), :]
m_vals_opti[i, :] = sol_state[7 + 7 * (i - 1), :]
end
sol_time_plot = LinRange(tvec[1], tvec[end], N)
x_moon_pos = [x_moon(t) for t in sol_time_plot]
y_moon_pos = [y_moon(t) for t in sol_time_plot]
z_moon_pos = [z_moon(t) for t in sol_time_plot]
x_earth_pos = [x_earth(t) for t in sol_time_plot]
y_earth_pos = [y_earth(t) for t in sol_time_plot]
z_earth_pos = [z_earth(t) for t in sol_time_plot]
p = plot(x_moon_pos, y_moon_pos, z_moon_pos, label = "moon orbit")
plot!(x_earth_pos, y_earth_pos, z_earth_pos, label = "earth orbit")
scatter!([x_vals_opti[1, 1]], [y_vals_opti[1, 1]], [z_vals_opti[1, 1]], label = "previous occultation")
for i in range(1, nb_arcs)
rand_color = RGB(rand(), rand(), rand())
plot!(x_vals_opti[i, :], y_vals_opti[i, :], z_vals_opti[i, :], size=(600, 600), camera = (30,30), xlabel="x", ylabel="y", zlabel="z", label="solar sail trajectory", linewidth = 2, color = "black", legend =:topleft)
scatter!([x_vals_opti[i, end]], [y_vals_opti[i, end]], [z_vals_opti[i, end]], label = "next occultation")
end
xlims!(-2, 2)
ylims!(-2, 2)
zlims!(-1.5, 1.5)
display(p)
# arrow_factor = 1e-1
# quiver!(x_vals_opti[u_norm .> 0.1][1:2:end], y_vals_opti[u_norm .> 0.1][1:2:end], z_vals_opti[u_norm .> 0.1][1:2:end], quiver=(ux[u_norm .> 0.1][1:2:end] * arrow_factor, uy[u_norm .> 0.1][1:2:end] * arrow_factor, uz[u_norm .> 0.1][1:2:end] * arrow_factor), label="thrust directions", arrow=true, color = "red")
end
This seems to work:
using SPICE
using OptimalControl
using NLPModelsIpopt
using DataInterpolations
using Downloads: download
using Plots
using Plots.PlotMeasures # for leftmargin, bottommargin
#using CSV
#using JLD2
# SPICE files that you have to download but only once
const LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
const SPK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp"
const PCK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc"
# Download kernels
download(LSK, "naif0012.tls")
download(SPK, "de440.bsp")
download(PCK, "pck00010.tpc")
furnsh("naif0012.tls") # Leap seconds kernel
furnsh("de440.bsp") # Ephemeris kernel
furnsh("pck00010.tpc") # PCK kernel for planetary constants
# Astrodynamical parameters
G = 6.6743e-11 / 1000^3 # km^3 kg^-1 s^-2
M_moon = 7.346e22 # kg
M_earth = 5.972168e24 # kg
M_sun = 1.989e30 # kg
AU = 1.495978707e11 / 1e3 # km
LU = 384399 # km
TU = sqrt(LU^3 / G / (M_moon + M_earth))
MUnit = M_moon + M_earth
mu_Earth = 3.9860044188e14 / 1000^3 / LU^3 * TU^2
mu_Moon = 4.90486959e12 / 1000^3 / LU^3 * TU^2
mu_Sun = 1.327124400189e20 / 1000^3 / LU^3 * TU^2
mass_moon = M_moon / MUnit
mass_earth = M_earth / MUnit
mass_sun = M_sun / MUnit
# Spacecraft parameters
g0 = 9.80665 / 1000 / LU * TU^2
Isp = 1660 / TU
wet_mass = 344 / MUnit
fuel_mass = 90 / MUnit
dry_mass = (wet_mass - fuel_mass)
wet_by_dry = (dry_mass + fuel_mass) / dry_mass
Tmax = 0.090 / 1000 / MUnit / LU * TU^2 #thrust force in N
Q = Tmax / Isp / g0
mass_scaling_factor = wet_mass
# Mission data that I just copy paste here for simplicity
t0_J200 = 1.107747244e9
tf_J200 = 1.112803048e9
N = 15000
t_values = LinRange(t0_J200, tf_J200, N-3)
t_values = vcat(-5e9, 0.0, collect(t_values), tf_J200 * 5)
states_moon = zeros(6, N)
states_earth = zeros(6, N)
states_sun = zeros(6, N)
t0 = t0_J200 / TU
tf = tf_J200 / TU
# Interpolations of the bodies positions/velocities
for i = 1:N
states_moon[:,i], _ = spkezr("MOON", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
states_earth[:,i], _ = spkezr("EARTH", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
states_sun[:,i], _ = spkezr("SUN", t_values[i], "ECLIPJ2000", "NONE", "EARTH MOON BARYCENTER")
end
# Interpolation of the position and velocity of the celestial bodies (Sun, Moon, Earth)
x_moon = DataInterpolations.LinearInterpolation(states_moon[1,:] / LU, t_values / TU)
y_moon = DataInterpolations.LinearInterpolation(states_moon[2,:] / LU, t_values / TU)
z_moon = DataInterpolations.LinearInterpolation(states_moon[3,:] / LU, t_values / TU)
x_earth = DataInterpolations.LinearInterpolation(states_earth[1,:] / LU, t_values / TU)
y_earth = DataInterpolations.LinearInterpolation(states_earth[2,:] / LU, t_values / TU)
z_earth = DataInterpolations.LinearInterpolation(states_earth[3,:] / LU, t_values / TU)
x_sun = DataInterpolations.LinearInterpolation(states_sun[1,:] / LU, t_values / TU)
y_sun = DataInterpolations.LinearInterpolation(states_sun[2,:] / LU, t_values / TU)
z_sun = DataInterpolations.LinearInterpolation(states_sun[3,:] / LU, t_values / TU)
# Times
t0 = 2952.5065962806448
t1 = 2953.955962786101
t2 = 2954.4423792964794
t3 = 2955.398985817685
t4 = 2957.4519685112937
t5 = 2959.0986010622196
t6 = 2960.6895789585565
# States at the previous times
ap0 = [-0.992463207586249; 0.7937735622182638; 0.07902677416947182; 0.4497665242298916; 0.5768992856452368; -0.08990077540870026]
occult1 = [0.024011456548615025; 1.0839463211562548; -0.06721110889958525; 0.8370532723530986; -0.2849379073937811; -0.09055054145691602]
lunar2 = [0.42191676845234877; 0.8365166413216278; -0.10058044976319964; 1.0641134073373073; -0.7420861511380829; 0.2221102307286784]
ap3 = [0.8866648294530993; 0.07703956973358668; 0.28758014926093095; 0.010909379448810359; -1.0104415437867642; 0.3065975859809964]
lunar4 = [-0.5240519542473246; -0.8727170323241018; 0.06983526263136841; -1.0656424593386877; 0.4013541226489389; -0.24177176981236737]
ap5 = [-1.2845625006191574; 0.07452729181167739; -0.02061966816776342; 0.03379116945381267; 0.7236740339189934; -0.05964259399261467]
occult6 = [-0.5063291057829847; 0.9375405520827109; -0.0856392904912281; 0.906825723975573; 0.17758145410610057; -0.009952663265456345]
# Changement de variable pour le temps
function phi(s, t0, t1)
return t0 + (t1 - t0) * s
end
# Dynamique
function Fu(t, r, v, m, u)
# Earth & Moon & Sun gravity
m = m * mass_scaling_factor
x_e = [x_earth(t), y_earth(t), z_earth(t)]
x_m = [x_moon(t), y_moon(t), z_moon(t)]
x_s = [x_sun(t), y_sun(t), z_sun(t)]
normrE = sqrt(sum((x_e - r).^2))
normrM = sqrt(sum((x_m - r).^2))
normrS = sqrt(sum((x_s - r).^2))
normrSE = sqrt(sum((x_e - x_s).^2))
normrMS = sqrt(sum((x_m - x_s).^2))
dv = mu_Earth / normrE^3 * (x_e - r) +
mu_Moon / normrM^3 * (x_m - r) +
mu_Sun / normrS^3 * (x_s - r) -
(1 / (mass_earth + mass_moon)) * (
x_s * (mu_Sun * mass_earth / normrSE^3 + mu_Moon * mass_sun / normrMS^3) -
x_m * (mu_Moon * mass_sun / normrMS^3) -
x_e * (mu_Sun * mass_earth / normrSE^3)
)
dvu = Tmax * u / m
dv = dv + dvu
dm = - Q * sqrt(sum(u.^2)) / mass_scaling_factor
dx = [v..., dv..., dm]
return dx
end
# Problem with 1 arc
@def ocp1arc begin
s ∈ [0, 1], time
y ∈ R^7, state
w ∈ R^3, control
r1 = y[1:3]
v1 = y[4:6]
m1 = y[7]
u1 = w[1:3]
[-5, -5, -5] ≤ r1(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v1(s) ≤ [5, 5, 5]
0.5 ≤ m1(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1, eq_u1
r1(0) == occult1[1:3]
v1(0) == occult1[4:6]
m1(0) == 1
r1(1) == lunar2[1:3]
v1(1) == lunar2[4:6]
ẏ(s) == (t2 - t1) * Fu(phi(s, t1, t2), r1(s), v1(s), m1(s), u1(s))
m1(1) → max
end
sol = solve(ocp1arc; grid_size = 100, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5)
sol = solve(ocp1arc; grid_size = 500, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, init = sol)
sol = solve(ocp1arc; grid_size = 900, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, init = sol)
plot(sol; state_bounds_style=:none, leftmargin=10mm, bottommargin=5mm)
# Problem with 3 arcs
@def ocp3arc begin
s ∈ [ 0, 1 ], time
y ∈ R^21, state
w ∈ R^9, control
# arc 1
x1 = y[1:7]
r1 = x1[1:3]
v1 = x1[4:6]
m1 = x1[7]
u1 = w[1:3]
# arc 2
x2 = y[8:14]
r2 = x2[1:3]
v2 = x2[4:6]
m2 = x2[7]
u2 = w[4:6]
# arc 3
x3 = y[15:21]
r3 = x3[1:3]
v3 = x3[4:6]
m3 = x3[7]
u3 = w[7:9]
# path constraints: arc 1
[-5, -5, -5] ≤ r1(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v1(s) ≤ [5, 5, 5]
0.5 ≤ m1(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1, eq_u1
# path constraints: arc 2
[-5, -5, -5] ≤ r2(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v2(s) ≤ [5, 5, 5]
0.5 ≤ m2(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u2(s).^2) ≤ 1, eq_u2
# path constraints: arc 3
[-5, -5, -5] ≤ r3(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v3(s) ≤ [5, 5, 5]
0.5 ≤ m3(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u3(s).^2) ≤ 1, eq_u3
# initial conditions: arc 1
r1(0) == ap0[1:3]
v1(0) == ap0[4:6]
m1(0) == 1
# final conditions: arc 1
r1(1) == occult1[1:3]
v1(1) == occult1[4:6]
# matching conditions: arcs 1 and 2
x2(0) - x1(1) == [0, 0, 0, 0, 0, 0, 0]
# final conditions: arc 2
r2(1) == lunar2[1:3]
v2(1) == lunar2[4:6]
# matching conditions: arcs 2 and 3
x3(0) - x2(1) == [0, 0, 0, 0, 0, 0, 0]
# final conditions: arc 3
r3(1) == ap3[1:3]
v3(1) == ap3[4:6]
# dynamics
ẏ(s) == [
(t1 - t0) * Fu(phi(s, t0, t1), r1(s), v1(s), m1(s), u1(s))...,
(t2 - t1) * Fu(phi(s, t1, t2), r2(s), v2(s), m2(s), u2(s))...,
(t3 - t2) * Fu(phi(s, t2, t3), r3(s), v3(s), m3(s), u3(s))...,
]
# objective
m3(1) → max
end
sol = solve(ocp3arc; grid_size = 100, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, disc_method=:gauss_legendre_2)
#sol = solve(ocp3arc; grid_size = 500, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, disc_method=:gauss_legendre_2, init = sol)
#sol = solve(ocp3arc; grid_size = 900, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, init = sol)
plot(sol; state_bounds_style=:none, leftmargin=10mm, bottommargin=5mm)
# Problem with 4 arcs
@def ocp4arc begin
s ∈ [ 0, 1 ], time
y ∈ R^28, state
w ∈ R^12, control
# arc 1
x1 = y[1:7]
r1 = x1[1:3]
v1 = x1[4:6]
m1 = x1[7]
u1 = w[1:3]
# arc 2
x2 = y[8:14]
r2 = x2[1:3]
v2 = x2[4:6]
m2 = x2[7]
u2 = w[4:6]
# arc 3
x3 = y[15:21]
r3 = x3[1:3]
v3 = x3[4:6]
m3 = x3[7]
u3 = w[7:9]
# arc 4
x4 = y[22:28]
r4 = x4[1:3]
v4 = x4[4:6]
m4 = x4[7]
u4 = w[10:12]
# path constraints: arc 1
[-5, -5, -5] ≤ r1(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v1(s) ≤ [5, 5, 5]
0.5 ≤ m1(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u1(s).^2) ≤ 1, eq_u1
# path constraints: arc 2
[-5, -5, -5] ≤ r2(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v2(s) ≤ [5, 5, 5]
0.5 ≤ m2(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u2(s).^2) ≤ 1, eq_u2
# path constraints: arc 3
[-5, -5, -5] ≤ r3(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v3(s) ≤ [5, 5, 5]
0.5 ≤ m3(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u3(s).^2) ≤ 1, eq_u3
# path constraints: arc 3
[-5, -5, -5] ≤ r4(s) ≤ [5, 5, 5]
[-5, -5, -5] ≤ v4(s) ≤ [5, 5, 5]
0.5 ≤ m4(s) ≤ 1 # careful here ! dry mass > 0
0 ≤ sum(u4(s).^2) ≤ 1, eq_u4
# initial conditions: arc 1
r1(0) == ap0[1:3]
v1(0) == ap0[4:6]
m1(0) == 1
# final conditions: arc 1
r1(1) == occult1[1:3]
v1(1) == occult1[4:6]
# matching conditions: arcs 1 and 2
x2(0) - x1(1) == [0, 0, 0, 0, 0, 0, 0]
# final conditions: arc 2
r2(1) == lunar2[1:3]
v2(1) == lunar2[4:6]
# matching conditions: arcs 2 and 3
x3(0) - x2(1) == [0, 0, 0, 0, 0, 0, 0]
# final conditions: arc 3
r3(1) == ap3[1:3]
v3(1) == ap3[4:6]
# matching conditions: arcs 3 and 4
x4(0) - x3(1) == [0, 0, 0, 0, 0, 0, 0]
# final conditions: arc 4
r4(1) == lunar4[1:3]
v4(1) == lunar4[4:6]
# dynamics
ẏ(s) == [
(t1 - t0) * Fu(phi(s, t0, t1), r1(s), v1(s), m1(s), u1(s))...,
(t2 - t1) * Fu(phi(s, t1, t2), r2(s), v2(s), m2(s), u2(s))...,
(t3 - t2) * Fu(phi(s, t2, t3), r3(s), v3(s), m3(s), u3(s))...,
(t4 - t3) * Fu(phi(s, t3, t4), r4(s), v4(s), m4(s), u4(s))...,
]
# objective
m4(1) → max
end
sol = solve(ocp4arc; grid_size = 100, max_iter = 100, tol = 1e-6, acceptable_tol = 1e-5, disc_method=:gauss_legendre_2)
#sol = solve(ocp4arc; grid_size = 500, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, disc_method=:gauss_legendre_2, init = sol)
#sol = solve(ocp4arc; grid_size = 900, max_iter = 3000, tol = 1e-6, acceptable_tol = 1e-5, init = sol)
plot(sol; state_bounds_style=:none, leftmargin=10mm, bottommargin=5mm)
Okay, donc le gros changement est de mettre les "..." dans les vecteurs. Merci beaucoup, je vais tester tout ça
Merci beaucoup @ocots ! It seems to be working indeed with the dots. Such a relief 😁
