[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