SimpleTrustRegion() diverges while SimpleNewtonRaphson() converges in MWE

[GrantHecht]

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()

[ChrisRackauckas]

@avik-pal @FHoltorf did SimpleTrustRegion get all of the updates that TrustRegion did?

[avik-pal]

I dont think so, the step size here is still

SimpleNonlinearSolve.jl/src/utils.jl

Line 63 in a65354a

δN = -H \ g

instead of $J$ \ $\delta f$

[FHoltorf]

Yes, it doesn't seem to have been fixed. The key issue seems to lie in

SimpleNonlinearSolve.jl/src/trustRegion.jl

Line 137 in a65354a

H = ∇f * ∇f

and

https://github.com/SciML/SimpleNonlinearSolve.jl/blob/a65354ab2617f6f9e7bbf9c30c712e7c8a1c97fc/src/trustRegion.jl#L191C24-L191C24

which incorrectly compute what $H$ and $g$ should be, i.e., $H = J'J$ and $g = J'f$. With those definitions, $-H \backslash g$ is technically correct (it just tends to be numerically unfavorable compared to computing $-J \backslash f$ if $J$ is full rank). I can go over it tomorrow and open a PR.

[ChrisRackauckas]

Fix in latest patch