From 12dbf5df3eae59752591006636e267810583597b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?= Date: Tue, 2 Jan 2024 07:29:37 -0600 Subject: [PATCH] Fixes to TSA --- .../orbit_determination.jl | 2 +- src/orbit_determination/tooshortarc.jl | 179 ++++++++++-------- src/propagation/parameters.jl | 8 +- 3 files changed, 102 insertions(+), 87 deletions(-) diff --git a/src/orbit_determination/orbit_determination.jl b/src/orbit_determination/orbit_determination.jl index 8665f585..bc2f66e9 100644 --- a/src/orbit_determination/orbit_determination.jl +++ b/src/orbit_determination/orbit_determination.jl @@ -44,7 +44,7 @@ function adaptative_maxsteps(radec::Vector{RadecMPC{T}}) where {T <: AbstractFlo return ceil(Int, (Δ_day + 360)/13) end =# - return 100 + return 500 end @doc raw""" diff --git a/src/orbit_determination/tooshortarc.jl b/src/orbit_determination/tooshortarc.jl index 0b58c538..c7c238f8 100644 --- a/src/orbit_determination/tooshortarc.jl +++ b/src/orbit_determination/tooshortarc.jl @@ -292,6 +292,9 @@ function propres(radec::Vector{RadecMPC{T}}, jd0::T, q0::Vector{U}, # Backward (forward) integration bwd = propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_bwd, q0, params) fwd = propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_fwd, q0, params) + if bwd.t[end] > t0 - jd0 || fwd.t[end] < tf - jd0 + return bwd, fwd, Vector{OpticalResidual{T, U}}(undef, 0) + end # Sun (Earth) ephemeris eph_su = selecteph(sseph, su) eph_ea = selecteph(sseph, ea) @@ -306,118 +309,128 @@ end @doc raw""" adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, - params::Parameters{T}; α::T = 25.0, β_1::T = 0.5, β_2::T = 0.85, - ϵ::T = 1e-8, Qtol::T = 0.01) where {T <: AbstractFloat} + params::Parameters{T}; η::T = 25.0, μ::T = 0.75, ν::T = 0.9, + ϵ::T = 1e-8, Qtol::T = 0.001) where {T <: AbstractFloat} + +Adaptative moment estimation (ADAM) minimizer of normalized mean square +residual over and admissible region `A`. -ADAM minimizer of root mean square error over `A`. +!!! warning + This function will set the (global) `TaylorSeries` variables to `δρ δvρ`. + +!!! reference + See Algorithm 1 of https://doi.org/10.48550/arXiv.1412.6980. """ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, - params::Parameters{T}; α::T = 25.0, β_1::T = 0.5, β_2::T = 0.85, - ϵ::T = 1e-8, Qtol::T = 0.01) where {T <: AbstractFloat} - # Origin - x0 = zeros(T, 2) + params::Parameters{T}; η::T = 25.0, μ::T = 0.75, ν::T = 0.9, + ϵ::T = 1e-8, Qtol::T = 0.001) where {T <: AbstractFloat} + # Initial time of integration [julian days] + jd0 = datetime2julian(A.date) # Scaling factors scalings = [A.ρ_domain[2] - A.ρ_domain[1], A.v_ρ_domain[2] - A.v_ρ_domain[1]] / 1_000 # Jet transport variables dq = scaled_variables("dρ dvρ", scalings, order = 1) + # Maximum number of iterations + maxiter = params.maxiter # Allocate memory - ρs = Vector{T}(undef, params.maxiter+1) - v_ρs = Vector{T}(undef, params.maxiter+1) - Qs = fill(T(Inf), params.maxiter+1) - # Initial time of integration [julian days] - jd0 = datetime2julian(A.date) - # Initial conditions - ρs[1] = ρ - v_ρs[1] = v_ρ - q0 = topo2bary(A, ρ + dq[1], v_ρ + dq[2]) - _, _, res = propres(radec, jd0, q0, params) - Q = nms(res) - Qs[1] = Q(x0) + ρs = Vector{T}(undef, maxiter+1) + v_ρs = Vector{T}(undef, maxiter+1) + Qs = fill(T(Inf), maxiter+1) + # Origin + x0 = zeros(T, 2) + # First momentum m = zeros(T, 2) - v = zeros(T, 2) - # Number of iterations - niter = 0 + _m_ = zeros(T, 2) + # Second momentum + n = zeros(T, 2) + _n_ = zeros(T, 2) # Gradient descent - while niter < params.maxiter - # Update number of iterations - niter += 1 + for t in 1:maxiter+1 + # Current position in admissible region + ρs[t] = ρ + v_ρs[t] = v_ρ + # Current barycentric state vector + q = topo2bary(A, ρ + dq[1], v_ρ + dq[2]) + # Propagation and residuals + _, _, res = propres(radec, jd0, q, params) + iszero(length(res)) && break + # Current Q + Q = nms(res) + Qs[t] = Q(x0) + # Convergence condition + t > 1 && abs(Qs[t] - Qs[t-1]) / Qs[t] < Qtol && break # Gradient of objective function - dQ = TaylorSeries.gradient(Q)(x0) - # Momentum - m = β_1 * m + (1 - β_1) * dQ - _m_ = m / (1 - β_1^niter) - # Sum of square of past gradients - v = β_2 * v + (1 - β_2) * dQ .^ 2 - _v_ = v / (1 - β_2^niter) + g_t = TaylorSeries.gradient(Q)(x0) + # First momentum + m .= μ * m + (1 - μ) * g_t + _m_ .= m / (1 - μ^t) + # Second momentum + n .= ν * n + (1 - ν) * g_t .^ 2 + _n_ .= n / (1 - ν^t) # Step - x1 = x0 - α * _m_ ./ (sqrt.(_v_) .+ ϵ) + x1 = x0 - η * _m_ ./ (sqrt.(_n_) .+ ϵ) # Update values - ρ, v_ρ = bary2topo(A, q0(x1)) + ρ, v_ρ = bary2topo(A, q(x1)) # Projection ρ, v_ρ = boundary_projection(A, ρ, v_ρ) - # New initial conditions - q0 = topo2bary(A, ρ + dq[1], v_ρ + dq[2]) - # Update obbjective function - _, _, res = propres(radec, jd0, q0, params) - Q = nms(res) - # Convergence condition - if abs(Q(x0) - Qs[niter]) < Qtol - break - end - # Save current variables - ρs[niter+1] = ρ - v_ρs[niter+1] = v_ρ - Qs[niter+1] = Q(x0) end + # Find point with smallest Q + t = argmin(Qs) + # Return path + return view(ρs, 1:t), view(v_ρs, 1:t), view(Qs, 1:t) +end - niter = argmin(Qs) +@doc raw""" + tsals(radec::Vector{RadecMPC{T}}, nights::Vector{ObservationNight{T}}, + jd0::T, q0::Vector{T}, params::Parameters{T}; maxiter::Int = 5, + Qtol::T = 0.1) where {T <: AbstractFloat} - return ρs[1:niter], v_ρs[1:niter], Qs[1:niter] -end +Used within [`tooshortarc`](@ref) to compute an orbit from a point in an +admissible region via least squares. -# Special method of tryls for tooshortarc -function tryls(radec::Vector{RadecMPC{T}}, nights::Vector{ObservationNight{T}}, - jd0::T, q0::Vector{T}, params::Parameters{T}; maxiter::Int = 5) where {T <: AbstractFloat} - # Allocate memory - sols = zeros(NEOSolution{T, T}, maxiter) - # Origin - x0 = zeros(T, 6) +!!! warning + This function will set the (global) `TaylorSeries` variables to `δx₁ δx₂ δx₃ δx₄ δx₅ δx₆`. +""" +function tsals(radec::Vector{RadecMPC{T}}, nights::Vector{ObservationNight{T}}, + jd0::T, q0::Vector{T}, params::Parameters{T}; maxiter::Int = 5, + Qtol::T = 0.1) where {T <: AbstractFloat} # Scaling factors scalings = abs.(q0) ./ 10^5 # Jet transport variables dq = scaled_variables("dx", scalings, order = 6) - # Number of iterations - niter = 1 - # Minimization loop - while niter <= maxiter + # Allocate memory + sols = zeros(NEOSolution{T, T}, maxiter+1) + Qs = fill(T(Inf), maxiter+1) + # Origin + x0 = zeros(T, 6) + # Least squares + for t in 1:maxiter+1 # Initial conditions q = q0 + dq # Propagation & residuals bwd, fwd, res = propres(radec, jd0, q, params) + iszero(length(res)) && break # Least squares fit - fit = tryls(res, x0, 5) + fit = tryls(res, x0, params.niter) # Case: unsuccessful fit if !fit.success || any(diag(fit.Γ) .< 0) break end # Current solution - sols[niter] = evalfit(NEOSolution(nights, bwd, fwd, res, fit, scalings)) - # Convergence condition - if niter > 1 && nrms(sols[niter-1]) - nrms(sols[niter]) < 0.1 - break + sols[t] = evalfit(NEOSolution(nights, bwd, fwd, res, fit, scalings)) + Qs[t] = nrms(sols[t]) + # Convergence conditions + if t > 1 + Qs[t] > Qs[t-1] && break + abs(Qs[t] - Qs[t-1]) < Qtol && break end # Update values q0 = q(fit.x) - # Update scaling factors - scalings = abs.(q0) ./ 10^5 - # Update Jet transport variables - dq = scaled_variables("dx", scalings, order = 6) - # Update number of iterations - niter += 1 end - - return sols[niter] - + # Find solution with smallest Q + t = argmin(Qs) + # Return solution + return sols[t] end # Order in which to check nights in tooshortarc @@ -449,7 +462,7 @@ function tooshortarc(radec::Vector{RadecMPC{T}}, nights::Vector{ObservationNight params::Parameters{T}) where {T <: AbstractFloat} # Allocate memory for output - sol = zero(NEOSolution{T, T}) + best_sol = zero(NEOSolution{T, T}) # Sort nights by tsanightorder idxs = sortperm(nights, lt = tsanightorder) @@ -461,16 +474,18 @@ function tooshortarc(radec::Vector{RadecMPC{T}}, nights::Vector{ObservationNight ρ = sum(A.ρ_domain) / 2 v_ρ = sum(A.v_ρ_domain) / 2 # Minimization over admissible region - ρs, v_ρs, Qs = adam(radec, A, ρ, v_ρ, params) + ρs, v_ρs, _ = adam(radec, A, ρ, v_ρ, params) # Barycentric initial conditions q0 = topo2bary(A, ρs[end], v_ρs[end]) # 6 variables least squares - sol = tryls(radec, nights, datetime2julian(A.date), q0, params; maxiter = 5) - - if nrms(sol) < 1.5 - return sol + sol = tsals(radec, nights, datetime2julian(A.date), q0, params; maxiter = 5) + # Update best solution + if nrms(sol) < nrms(best_sol) + best_sol = sol + # Break condition + nrms(sol) < 1.5 && break end end - return sol + return best_sol end \ No newline at end of file diff --git a/src/propagation/parameters.jl b/src/propagation/parameters.jl index 3df2a9d4..72dda533 100644 --- a/src/propagation/parameters.jl +++ b/src/propagation/parameters.jl @@ -32,7 +32,7 @@ Parameters for all orbit determination functions. # Propagation Parameters -- `maxsteps::Int`: maximum number of steps for the integration (default: `100`). +- `maxsteps::Int`: maximum number of steps for the integration (default: `500`). - `μ_ast::Vector{T}`: vector of gravitational parameters (default: `μ_ast343_DE430[1:end]`). - `order::Int`: order of Taylor expansions wrt time (default: 25). - `abstol::T`: absolute tolerance used to compute propagation timestep (default: `1e-20`). @@ -57,16 +57,16 @@ Parameters for all orbit determination functions. # Admissible Region Parameters -- `maxiter::Int`: maximum number of iterations for admissible region `ADAM` optimizer (default: `200`). +- `maxiter::Int`: maximum number of iterations for admissible region `ADAM` optimizer (default: `100`). # Outlier Rejection Parameters - `max_per::T`: maximum allowed rejection percentage (default: `18.0`). """ -function Parameters(; maxsteps::Int = 100, μ_ast::Vector{T} = μ_ast343_DE430[1:end], +function Parameters(; maxsteps::Int = 500, μ_ast::Vector{T} = μ_ast343_DE430[1:end], order::Int = 25, abstol::T = 1e-20, parse_eqs::Bool = true, debias_table::String = "2018", niter::Int = 5, max_triplets::Int = 10, - varorder::Int = 5, Q_max::T = 5.0, maxiter::Int = 200, + varorder::Int = 5, Q_max::T = 5.0, maxiter::Int = 100, max_per::T = 18.0) where {T <: AbstractFloat} # Unfold debiasing matrix mpc_catalogue_codes_201X, truth, resol, bias_matrix = select_debiasing_table(debias_table)