First, I want to say thank you for developing this package! I've been using NLsolve.jl for years now and am excited by the increased performance I'll obtain with your fast implementations of NewtonRaphson and TrustRegion algorithms.
I've put together an MWE that involves finding the solution to a spacecraft trajectory optimization problem using simple forward shooting. I already know a solution to this problem, and am setting an initial guess that is very near this solution. As expected, SimpleNewtonRaphson() successfully converges to the solution, but SimpleTrustRegion() diverges.
I suppose this could be a non-issue as there are no convergence guarantees when solving general nonlinear systems, but from my experience with this problem, I find it strange that the trust region algorithm fails to converge in this case. I understand that SimpleTrustRegion() was implemented recently, so I wanted to post this issue in the hopes that it will help identify any bug that could be causing this behavior.
using DifferentialEquations
using StaticArrays
using SimpleNonlinearSolve
# State/Costate equations of motion for min-fuel trajectory optimization in CRTBP model
function crtbp_min_fuel_eom(x, p, t)
# Grab states
rx = x[1]; ry = x[2]; rz = x[3]
vx = x[4]; vy = x[5]; vz = x[6]
m = x[7]
# Grab costates
lam_rx = x[8]; lam_ry = x[9]; lam_rz = x[10]
lam_vx = x[11]; lam_vy = x[12]; lam_vz = x[13]
# Grab parameters
mu = p[1]
T_max = p[2]
c = p[3]
u = p[4]
return @fastmath begin
# Start of MATLAB generated code
t2 = mu+rx;
t3 = lam_vx*lam_vx;
t4 = lam_vy*lam_vy;
t5 = lam_vz*lam_vz;
t8 = ry*ry;
t9 = rz*rz;
t10 = 1.0/m;
t11 = mu-1.0;
t12 = t2-1.0;
t13 = t2*t2;
t17 = t3+t4+t5;
t14 = t12*t12;
t18 = t8+t9+t13;
t20 = 1.0/sqrt(t17);
t19 = t8+t9+t14;
# We can make this code faster so we will...
#t21 = 1.0/pow(t18,3.0/2.0);
#t22 = 1.0/pow(t18,5.0/2.0);
tt1 = sqrt(t18)
tt13 = tt1*tt1*tt1
tt15 = tt13*tt1*tt1
t21 = 1.0/tt13
t22 = 1.0/tt15
# Same thing here...
#t23 = 1.0/pow(t19,3.0/2.0);
#t24 = 1.0/pow(t19,5.0/2.0);
tt2 = sqrt(t19)
tt23 = tt2*tt2*tt2
tt25 = tt23*tt2*tt2
t23 = 1.0/tt23
t24 = 1.0/tt25
t26 = t11*t21;
t29 = ry*rz*t11*t22*3.0;
t25 = mu*t23;
t28 = mu*ry*rz*t24*3.0;
t27 = -t25;
return SA[
vx, vy, vz,
rx+vy*2.0+t2*t26+t12*t27-T_max*lam_vx*t10*t20*u,
ry-vx*2.0+ry*t26+ry*t27-T_max*lam_vy*t10*t20*u,
rz*t26+rz*t27-T_max*lam_vz*t10*t20*u,
-(T_max*u)/c,
-lam_vy*(mu*ry*t12*t24*3.0-ry*t2*t11*t22*3.0)-lam_vz*(mu*rz*t12*t24*3.0-rz*t2*t11*t22*3.0)-lam_vx*(t26+t27+mu*t14*t24*3.0-t11*t13*t22*3.0+1.0),
-lam_vz*(t28-t29)-lam_vx*(mu*ry*t12*t24*3.0-ry*t2*t11*t22*3.0)-lam_vy*(t26+t27+mu*t8*t24*3.0-t8*t11*t22*3.0+1.0),
-lam_vy*(t28-t29)-lam_vx*(mu*rz*t12*t24*3.0-rz*t2*t11*t22*3.0)+lam_vz*(t25-t26-mu*t9*t24*3.0+t9*t11*t22*3.0),
-lam_rx+lam_vy*2.0,
-lam_ry-lam_vx*2.0,
-lam_rz,
-(T_max*(t10*t10)*u)/t20,
]
end
end
function condition(u,t,integ)
c = integ.p[3]
Ξ»v = sqrt(u[11]*u[11] + u[12]*u[12] + u[13]*u[13])
return -c*Ξ»v/u[7] - u[14] + 1.0
end
function affect!(integrator)
params_pre_affect = integrator.p
new_throttle_factor = params_pre_affect[end] == 1 ? 0 : 1
integrator.p = (params_pre_affect[1], params_pre_affect[2], params_pre_affect[3], new_throttle_factor)
return nothing
end
function initialize!(c, u, t, integrator)
c = integrator.p[3]
Ξ»v = sqrt(u[11]*u[11] + u[12]*u[12] + u[13]*u[13])
S = -c*Ξ»v/u[7] - u[14] + 1.0
sf = 1
if S > 0
sf = 0
end
ps = integrator.p
integrator.p = (ps[1], ps[2], ps[3], sf)
return nothing
end
function shooting_function(Ξ»0, p)
# Grab parameters
x0 = p[1]
xf = p[2]
tspan = p[3]
mu = p[4]
T_max = p[5]
c = p[6]
# Form full initial state
y0 = SA[
x0[1], x0[2], x0[3], x0[4], x0[5], x0[6], x0[7],
Ξ»0[1], Ξ»0[2], Ξ»0[3], Ξ»0[4], Ξ»0[5], Ξ»0[6], Ξ»0[7],
]
# Integrate dynamics from initial condition
cb = ContinuousCallback(
condition,
affect!;
initialize = initialize!,
rootfind = SciMLBase.RightRootFind,
)
prob = ODEProblem{false,SciMLBase.FullSpecialize}(
crtbp_min_fuel_eom,
y0,
tspan,
(mu, T_max, c, 1),
)
sol = solve(
prob,
Vern9();
abstol = 1e-14,
reltol = 1e-14,
save_everystep = false,
save_start = false,
initialize_save = false,
maxiters = 1e7,
callback = cb,
)
sxf = sol[end]
# Return final boundary condition residual
return SA[
sxf[1] - xf[1],
sxf[2] - xf[2],
sxf[3] - xf[3],
sxf[4] - xf[4],
sxf[5] - xf[5],
sxf[6] - xf[6],
sxf[14],
]
end
function main()
# CRTBP parameters
G = 6.673e-20
m_e = 5.9742e24 # kg
m_m = 7.3483e22 # kg
r_em = 384400.0 # km
mu = 1.21506038e-2
# Define spacecraft parameters
T_max_si = 10e-3 # kg*km/s^2
I_sp_si = 3000.0 # s
g_0_si = 9.81e-3 # km/s^2
c_si = I_sp_si * g_0_si
m0_si = 1500.0 # kg
# Units
LU = r_em
TU = 1.0 / sqrt((G*(m_e + m_m)) / LU^3)
MU = m0_si
# Scale parameters
T_max_nd = T_max_si * TU^2 / (MU * LU)
c_nd = c_si * TU / LU
# Define time of flight
tspan = (0.0, 8.6404*86400.0 / TU)
# Define initial spacecraft state
x0 = SA[
-0.0194885115, -0.0160334798, 0.0,
8.9188819237, -4.0817936888, 0.0, 1.0,
]
# Define initial co-state guess (Traj. B min. fuel solution from Acta Astronautica paper)
Ξ»0 = SA[
15.68499615832586,
33.031211178919726,
-0.09381690521953079,
-0.10207501688908052,
0.045012820414046514,
-0.00015063964363906,
0.1334329267459845,
]
# Define target spacecraft state (note final mass unconstrainted)
xf = SA[
0.8233851820, 0.0, -0.0222775563,
0.0, 0.1341841703, 0.0,
]
# Define Nonlinear Problem
params = (x0, xf, tspan, mu, T_max_nd, c_nd)
# Try to solve
prob = NonlinearProblem{false}(shooting_function, Ξ»0, params)
# Here, SimpleTrustRegion() diverges whereas SimpleNewtonRaphson() converges
Ξ»0 = solve(prob, SimpleTrustRegion(); verbose = true) # Unfortunately, verbose = true does nothing currently :(
#Ξ»0 = solve(prob, SimpleNewtonRaphson(); verbose = true)
# Compute residual
resid = shooting_function(Ξ»0, params)
# If residual is small enough, we have a solution
max_resid = maximum(abs.(resid))
display(max_resid)
return nothing
end
main()
![dependabot[bot] avatar](https://avatars.githubusercontent.com/in/29110?v=4 "dependabot[bot]")
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent