[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.

```julia
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
```