28 dev kepler add examples

.gitignore

@@ -30,3 +30,6 @@ Manifest.toml
 **/Manifest.toml
 # ...except the one inside docs/src/assets
 !docs/src/assets/Manifest.toml
+
+# Ignore references 
+references/

save/3body_L1.jl

@@ -0,0 +1,79 @@
+using OptimalControl
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+#using MINPACK
+using ForwardDiff
+using LinearAlgebra
+
+ε = 0.01 
+
+## Problem definition
+
+G = 6.61e10-11
+m₁ = 5.972e24          # Mass of the Earth
+m₂ = 7.348e22          # Mass of the Moon
+μ = m₂ / (m₁ + m₂)     # Ratio of primary masses
+M = 1000               # Mass of the satellite (kg)
+# !!!!!!!!!!!!
+γ_max = 1              # Poussée max
+T_max = 1              # Maximal thrust
+# !!!!!!!!!!!!
+
+# 1st body is located at (-μ, 0)
+# 2nd body is located at (1-μ, 0)
+
+# Kepler law speed in circular movement: v² = G * M / d
+d_earth = 0.1
+d_moon = 0.1
+
+v_earth_circular = sqrt(G * m₁ / d_earth)
+v_moon_circular = sqrt(G * m₂ / d_moon)
+
+q0 = [-μ - d_earth, 0.0]   # Initial position (km) (e.g. GEO Earth)
+v0 = [0, v_earth_circular]       # Initial velocity (km/s)
+qf = [1 - μ + d_moon, 0.0]    # Target position (km) (e.g. GEO Moon)
+vf = [0.0, -v_moon_circular]       # Target velocity (km/s)
+x0 = vcat(q0, v0)
+xf = vcat(qf, vf)
+
+tf_guess = 20000.0          # Time guess in seconds (5.5 hours)
+
+# =============== INITIAL GUESS ===============
+x(t) = x0 + (xf - x0) * t / tf_guess
+u(t) = [0.01, 0.01] # initial control guess
+
+
+nlp_init = (state = x, control = u, variable = tf_guess)
+
+# =============== OCP DEFINITION ===============
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x ∈ R⁴, state
+    u ∈ R², control
+    x(0) == x0
+    x(tf) == xf
+    r₁₃ = sqrt((x₁(t) + μ)^2 + x₂(t)^2) + 1e-9  # Avoid division by zero
+    r₂₃ = sqrt((x₁(t) - 1 + μ)^2 + x₂(t)^2) + 1e-9  # Avoid division by zero
+    F0 = [  x₃(t), 
+            x₄(t), 
+            2*x₄(t) + x₁(t) - ((1 - μ) / r₁₃^3)*(x₁(t) + μ) - (μ / r₂₃^3)*(x₁(t) - 1 + μ),
+            -2*x₃(t) + x₂(t) - ((1 - μ) / r₁₃^3)*x₂(t) - (μ / r₂₃^3)*x₂(t) ]
+    F1 = [0, 0, 1, 0]
+    F2 = [0, 0, 0, 1]
+    ẋ(t) == F0 + (T_max/M)*F1*u1(t) + (T_max/M)*F2*u₂(t)
+    u₁(t)^2 + u₂(t)^2 ≤ 1
+    # Fonction objectif avec régularisation pour éviter les problèmes de domaine
+    u_norm = sqrt(u₁(t)^2 + u₂(t)^2) + 1e-12  # Évite log(0)
+    u_complement = max(1e-12, 1 - u_norm)  # Évite log(≤0)
+    ∫(u_norm - (1-ε)*log(u_norm) - (1-ε)*log(u_complement)) → min
+end
+
+
+# =============== SOLVE ===============
+
+nlp_sol = OptimalControl.solve(ocp; init=nlp_init, grid_size=50, print_level=3)
+plot(nlp_sol, vars=:x, title="3-body Controlled Transfer")
+
+

save/three_body_time.jl

@@ -0,0 +1,53 @@
+using OptimalControl
+using CTBase
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+using ForwardDiff
+using LinearAlgebra
+
+## Problem definition (CR3BP)
+
+μ = 0.0121505856              # Earth-Moon reduced mass ratio
+ϵ = 1.0                       # Maximum thrust magnitude (normalized)
+q0 = [1.0 - μ - 0.01, 0.0]    # Initial position (near Earth)
+v0 = [0.0, 0.0]               # Initial velocity
+qf = [μ + 0.05, 0.0]          # Final position (near Moon)
+vf = [0.0, 0.0]               # Final velocity
+x0 = vcat(q0, v0)
+xf = vcat(qf, vf)
+tf_guess = 5.0               # Initial guess for final time
+
+# Gradient of the effective potential
+function grad_omega(q)
+    x, y = q
+    r1 = sqrt((x + μ)^2 + y^2) + 1e-6
+    r2 = sqrt((x - 1 + μ)^2 + y^2) + 1e-6
+    dΩdx = x - (1 - μ)*(x + μ)/r1^3 - μ*(x - 1 + μ)/r2^3
+    dΩdy = y - (1 - μ)*y/r1^3 - μ*y/r2^3
+    return [dΩdx, dΩdy]
+end
+
+
+x(t) = x0 + (xf - x0) * t / tf_guess
+u(t) = [0.01, 0.01] 
+nlp_init = (state = x, control = u, variable = tf_guess)
+
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    x = (q1, q2, v1, v2) ∈ R⁴, state
+    u ∈ R², control
+    x(0) == x0
+    x(tf) == xf
+    r1 = sqrt((q1(t) + μ)^2 + q2(t)^2) + 1e-6
+    r2 = sqrt((q1(t) - 1 + μ)^2 + q2(t)^2) + 1e-6
+    dΩ1 = q1(t) - (1 - μ)*(q1(t) + μ)/r1^3 - μ*(q1(t) - 1 + μ)/r2^3
+    dΩ2 = q2(t) - (1 - μ)*q2(t)/r1^3 - μ*q2(t)/r2^3
+    ẋ(t) == [v1(t), v2(t), 2v2(t) + dΩ1 + ϵ*u₁(t), -2v1(t) + dΩ2 + ϵ*u₂(t)]
+    u₁(t)^2 + u₂(t)^2 ≤ 1
+    tf → min
+end
+
+nlp_sol = OptimalControl.solve(ocp; init=nlp_init, grid_size=100, print_level=3)
+plot(nlp_sol, vars=:x, title="CR3BP - Minimum Time Transfer")

save/two_body.jl

@@ -0,0 +1,138 @@
+using OptimalControl
+using NLPModelsIpopt
+using OrdinaryDiffEq
+using Plots
+using MINPACK
+using ForwardDiff
+using LinearAlgebra
+
+Tmax = 60                                  # Maximum thrust in Newtons
+cTmax = 3600^2 / 1e6; T = Tmax * cTmax     # Conversion from Newtons to kg x Mm / h²
+mass0 = 1500                               # Initial mass of the spacecraft
+β = 1.42e-02                               # Engine specific impulsion
+μ = 5165.8620912                           # Earth gravitation constant
+P0 = 11.625                                # Initial semilatus rectum
+ex0, ey0 = 0.75, 0                         # Initial eccentricity
+hx0, hy0 = 6.12e-2, 0                      # Initial ascending node and inclination
+L0 = π                                     # Initial longitude
+Pf = 42.165                                # Final semilatus rectum
+exf, eyf = 0, 0                            # Final eccentricity
+hxf, hyf = 0, 0                            # Final ascending node and inclination
+ε = 1e-1                                   # Regularization parameter for logarithmic barrier
+
+asqrt(x; ε=1e-9) = sqrt(sqrt(x^2 + ε^2))
+function F0(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = zeros(eltype(x), 6)
+    F[6] = w^2 / (P * pdm)
+    return F
+end
+
+function F1(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    F = zeros(eltype(x), 6)
+    F[2] = pdm *   sl
+    F[3] = pdm * (-cl)
+    return F
+end
+
+function F2(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    F = zeros(eltype(x), 6)
+    F[1] = pdm * 2 * P / w
+    F[2] = pdm * (cl + (ex + cl) / w)
+    F[3] = pdm * (sl + (ey + sl) / w)
+    return F
+end
+
+function F3(x)
+    P, ex, ey, hx, hy, L = x
+    pdm = asqrt(P / μ)
+    cl = cos(L)
+    sl = sin(L)
+    w = 1 + ex * cl + ey * sl
+    pdmw = pdm / w
+    zz = hx * sl - hy * cl
+    uh = (1 + hx^2 + hy^2) / 2
+    F = zeros(eltype(x), 6)
+    F[2] = pdmw * (-zz * ey)
+    F[3] = pdmw *   zz * ex
+    F[4] = pdmw *   uh * cl
+    F[5] = pdmw *   uh * sl
+    F[6] = pdmw *   zz
+    return F
+end
+
+tf = 15
+Lf = 3π
+x0 = [P0, ex0, ey0, hx0, hy0, L0]
+xf = [Pf, exf, eyf, hxf, hyf, Lf]
+x(t) = x0 + (xf - x0) * t / tf
+u = [0.1, 0.5, 0.]
+nlp_init = (state=x, control=u, variable=tf)
+
+function min_tf()
+    @def ocp begin
+        tf ∈ R, variable
+        t ∈ [0, tf], time
+        x = (P, ex, ey, hx, hy, L) ∈ R⁶, state
+        u ∈ R³, control
+        x(0) == x0
+        x[1:5](tf) == xf[1:5]
+        mass = mass0 - β * T * t
+        ẋ(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
+        u₁(t)^2 + u₂(t)^2 + u₃(t)^2 ≤ 1
+        tf → min
+    end
+    return ocp
+end
+
+function min_conso(e)
+    @def ocp begin
+        t ∈ [0, tf], time
+        x = (P, ex, ey, hx, hy, L) ∈ R⁶, state
+        u ∈ R³, control
+        x(0) == x0
+        x[1:5](tf) == xf[1:5]
+        mass = mass0 - β * T * t
+        u_norm = sqrt(u₁(t)^2 + u₂(t)^2 + u₃(t)^2)
+        ẋ(t) == F0(x(t)) + T / mass * (u₁(t) * F1(x(t)) + u₂(t) * F2(x(t)) + u₃(t) * F3(x(t)))
+        1e-3 ≤ u_norm^2 ≤ 1
+        u_g = max(1e-10, 1 - u_norm)
+        ∫(u_norm - e * (log(u_norm) + log(u_g))) → min
+    end
+    return ocp
+end
+
+ocp1 = min_tf()
+sol = solve(ocp1; init=nlp_init, grid_size=100)
+Tmax = 100.0
+T = Tmax * cTmax
+tf = variable(sol)
+ocp = min_conso(ε)
+nlp_init = (state = state(sol), control = control(sol))
+
+eps = 1e-1
+i = nlp_init
+s = sol
+ocp = min_conso(eps)
+global s = solve(ocp; init=i, grid_size=500)
+global i = s
+
+
+plot(s; label="ε = $eps")
+
+t = time_grid(s)
+u = control(s)
+plot(t, norm∘u; label="‖u‖", xlabel="t")