TMJets crashing iteratively, what to try?

[willsharpless]

Describe the bug

I am trying to solve the reachable sets for an N-dimensional L96 model but unfortunately the TMJets algorithms causes the REPL to exit abruptly or the evaluations don't return (killed after 30 mins).

In the MWE below, N=9 but I have also tried with N=4 and see the same results. Note, the model includes disturbance but I have it set to 0 by default below as it wouldn't solve in either case. Moreover, after initial failure, I chose a very small X0 as an inf-ball of r = 0.01 with center defined by the natural evolution of the system after t=1 (solved by DifferentialEquations.jl) and a small tspan = (0., 0.1) for the reachability analysis to simplify the problem but to no avail.

I tried a series of orderQ's and orderT's seen in comments at the bottom of the MWE,

## Terminal Process is terminated with exit code 139:
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets()) # default (docs say) areo rderT=6, orderQ = 2
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets21a(abstol=1e-10))
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=4))

## doesn't crash but doesn't finish
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=8)) 
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=10, orderQ=12)) 
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-5, orderT=10, orderQ=12, adaptive=false)) 

starting from the default and then increasing, but none yielded results, although some caused the REPL to terminate and some just ran for a long time (30 min~) before I killed them.

Do you all have any ideas what might be causing the difficulty of the solution for this model? Do you have any thoughts on other things I might try? I also tried CARLIN but ran into a syntax bug or something.

L96 is well known to be chaotic but we are considering a very small timescale and space (w/o disturbance too), and for the most part the system trajectories should be near the limit cycle since I have let the system evolve for a second before this analysis.

Finally, I may say that adding @taylorize hasn't appeared to change anything but I understand it is sensitive to the fn structure.

To Reproduce

using LinearAlgebra, Plots
using ReachabilityAnalysis

N = 9; 
nd = N; 

function L96!(dx, x, p, t)
    F = 8.0 * ones(N)

    ## L96 states
    dx[1] = (x[2] - x[8]) * x[9] - x[1] + F[1] + x[N+1]
    dx[2] = (x[3] - x[9]) * x[1] - x[2] + F[2] + x[N+2]
    dx[3] = (x[4] - x[1]) * x[2] - x[3] + F[3] + x[N+3]
    dx[4] = (x[5] - x[2]) * x[3] - x[4] + F[4] + x[N+4]
    dx[5] = (x[6] - x[3]) * x[4] - x[5] + F[5] + x[N+5]
    dx[6] = (x[7] - x[4]) * x[5] - x[6] + F[6] + x[N+6]
    dx[7] = (x[8] - x[5]) * x[6] - x[7] + F[7] + x[N+7]
    dx[8] = (x[9] - x[6]) * x[7] - x[8] + F[8] + x[N+8]
    dx[9] = (x[1] - x[7]) * x[8] - x[9] + F[9] + x[N+9]
    
    ## Disturbance Dummy States
    dx[N+1:end] = zero(x[N+1:end])

    return dx
end

### Zonotoping 
t = 1.
Xeq = [6.84, 4.89, 4.13, 5.56, 6.51, 4.71, 3.06, 4.06, 6.16]

## Target
r𝒯 = 0.01;
c𝒯 = vcat(Xeq); Q𝒯 = inv(r𝒯) * diagm(ones(N))

X0 = Hyperrectangle(; low = c𝒯 - diag(inv(Q𝒯)), high = c𝒯 + diag(inv(Q𝒯)))
D0 = Hyperrectangle(; low = zeros(nd), high = zeros(nd))

X̂0auto = X0 × D0;
sys_auto = InitialValueProblem(BlackBoxContinuousSystem(L96!, length(low(X̂0auto))), X̂0auto);

# Terminal Process is terminated with exit code 139:
# sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets()) # default should (docs say)
# sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets21a(abstol=1e-10))
# sol = solve(sys, tspan=(t, 0.); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=4))

# doesn't crash but doesn't finish
sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=8)) 
# sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-10, orderT=10, orderQ=12)) 
# sol = solve(sys_auto, tspan=(0., 0.1); alg=TMJets(; maxsteps=1_000, abstol=1e-5, orderT=10, orderQ=12, adaptive=false)) 

[mforets]

If the system is 9-dimensional, don't you have out of bounds in x[N+1] and other terms?

EDIT: Ah, you pass X0 × D0 as initial states so that should be fine.

[mforets]

I'm getting some results with a slightly modified version:

Code:

using LinearAlgebra, Plots
using ReachabilityAnalysis

N = 9;
nd = N;

@taylorize function L96!(dx, x, p, t)
    # F = 8.0 * ones(N)

    ## L96 states
    dx[1] = (x[2] - x[8]) * x[9] - x[1] + 8.0 + x[N + 1]
    dx[2] = (x[3] - x[9]) * x[1] - x[2] + 8.0 + x[N + 2]
    dx[3] = (x[4] - x[1]) * x[2] - x[3] + 8.0 + x[N + 3]
    dx[4] = (x[5] - x[2]) * x[3] - x[4] + 8.0 + x[N + 4]
    dx[5] = (x[6] - x[3]) * x[4] - x[5] + 8.0 + x[N + 5]
    dx[6] = (x[7] - x[4]) * x[5] - x[6] + 8.0 + x[N + 6]
    dx[7] = (x[8] - x[5]) * x[6] - x[7] + 8.0 + x[N + 7]
    dx[8] = (x[9] - x[6]) * x[7] - x[8] + 8.0 + x[N + 8]
    dx[9] = (x[1] - x[7]) * x[8] - x[9] + 8.0 + x[N + 9]

    ## Disturbance Dummy States
    local zerox = zero(x[1])
    dx[10] = zerox
    dx[11] = zerox
    dx[12] = zerox
    dx[13] = zerox
    dx[14] = zerox
    dx[15] = zerox
    dx[16] = zerox
    dx[17] = zerox
    dx[18] = zerox

    return dx
end

### Zonotoping 
t = 1.0
Xeq = [6.84, 4.89, 4.13, 5.56, 6.51, 4.71, 3.06, 4.06, 6.16]

## Target
r𝒯 = 0.01;
c𝒯 = vcat(Xeq);
Q𝒯 = inv(r𝒯) * diagm(ones(N));

X0 = Hyperrectangle(; low=c𝒯 - diag(inv(Q𝒯)), high=c𝒯 + diag(inv(Q𝒯)))
D0 = Hyperrectangle(; low=zeros(nd), high=zeros(nd))

X̂0auto = X0 × D0;
sys_auto = InitialValueProblem(BlackBoxContinuousSystem(L96!, length(low(X̂0auto))), X̂0auto);

# Adjust abstol to tradeoff speed and accuracy.
@time sol = solve(sys_auto; tspan=(0.0, 1.0), alg=TMJets(; abstol=1e-14));
41.484481 seconds (153.99 M allocations: 55.377 GiB, 18.00% gc time)

julia> plot(sol, vars=(1, 2), xlab="x1", ylab="x2")

julia> plot(sol, vars=(0, 3), xlab="t", ylab="x3")

EDIT: to make the loop in the function generic, you can otherwise use a loop:

using LinearAlgebra, Plots
using ReachabilityAnalysis

const N = 9;
const nd = N;
const F = 8 * ones(N)

@taylorize function L96!(dx, x, p, t)
    ## L96 states
    dx[1] = (x[2] - x[8]) * x[9] - x[1] + F[1] + x[N + 1]
    dx[2] = (x[3] - x[9]) * x[1] - x[2] + F[2] + x[N + 2]
    dx[3] = (x[4] - x[1]) * x[2] - x[3] + F[3] + x[N + 3]
    dx[4] = (x[5] - x[2]) * x[3] - x[4] + F[4] + x[N + 4]
    dx[5] = (x[6] - x[3]) * x[4] - x[5] + F[5] + x[N + 5]
    dx[6] = (x[7] - x[4]) * x[5] - x[6] + F[6] + x[N + 6]
    dx[7] = (x[8] - x[5]) * x[6] - x[7] + F[7] + x[N + 7]
    dx[8] = (x[9] - x[6]) * x[7] - x[8] + F[8] + x[N + 8]
    dx[9] = (x[1] - x[7]) * x[8] - x[9] + F[9] + x[N + 9]

    ## Disturbance Dummy States
    local zerox = zero(x[1])
    for i in (N + 1):(2N)
        dx[i] = zerox
    end

    return dx
end

### Zonotoping 
t = 1.0
Xeq = [6.84, 4.89, 4.13, 5.56, 6.51, 4.71, 3.06, 4.06, 6.16]

## Target
r𝒯 = 0.01;
c𝒯 = vcat(Xeq);
Q𝒯 = inv(r𝒯) * diagm(ones(N));

X0 = Hyperrectangle(; low=c𝒯 - diag(inv(Q𝒯)), high=c𝒯 + diag(inv(Q𝒯)))
D0 = Hyperrectangle(; low=zeros(nd), high=zeros(nd))

X̂0auto = X0 × D0;
sys_auto = InitialValueProblem(BlackBoxContinuousSystem(L96!, length(low(X̂0auto))), X̂0auto);

# Terminal Process is terminated with exit code 139:
@time sol = solve(sys_auto; tspan=(0.0, 1.0), alg=TMJets(; abstol=1e-12));

plot(sol; vars=(1, 2))

[willsharpless]

Hey wow looks great! Works for me too. Seems like the dx[N+1:end] = zero(x[N+1:end]) call is what broke it. Apologies for not considering this, thank you!

[schillic]

I reopen here because this was tagged as a "bug". @mforets, is there some actionable point or should we change the tag?