Warning when trying to solve nonlinear ODE

[SebastianGuadalupe]

I tried to solve the warning changing max_steps and abs_tol but it didn't work

using ReachabilityAnalysis
@taylorize function cartpole!(du, u, p, t)
    local f, m, l, mt, g = 10, 0.1, 0.5, 1.1, 9.8
    
    du[1] = u[2]
    du[2] = (f + m*l*u[4]^2*sin(u[3])) / mt - ((m*l*(g*sin(u[3])- cos(u[3]))*((f+m*l*u[4]^2*sin(u[3]))/mt)) / l*(4/3 - m*cos(u[3])^2/mt)) * cos(u[3])/mt
    du[3] = u[4]
    du[4] = (g*sin(u[3])-cos(u[3])*((f+m*l*u[4]^2*sin(u[3]))/mt)) / l*(4/3 - m*cos(u[3])^2/mt) * cos(u[3])/mt
    
    return du
end

# define the initial-value problem
X₀ = Hyperrectangle(low=[-0.05, -0.05, -0.05, -0.05], high=[0.05, 0.05, 0.05, 0.05])

prob = @ivp(x' = cartpole!(x), dim: 4, x(0) ∈ X₀)

# solve it
sol = solve(prob, T=12.0, max_steps=200000, abs_tol=1e-20);

┌ Warning: Maximum number of integration steps reached; exiting.
└ @ ReachabilityAnalysis /home/sguadalupe/.julia/dev/ReachabilityAnalysis/src/Algorithms/TMJets/validated_integ.jl:517

[mforets]

The options max_steps and abs_tol are not passed to the solver; you have to use TMJets(max_steps=200000, abs_tol=1e-20) (that max_steps it too large though; we should see if there is another way eg. by splitting the initial conditions). That said, i'm considering to pass the kwargs... to the solver to avoid this kind of problem.

[mforets]

BTW note that the solution structure contains the options that were used, eg.

julia> sol = solve(prob, T=1.0, alg=TMJets(max_steps=10_000, abs_tol=1e-10));

julia> sol.alg
TMJets{Float64,ReachabilityAnalysis.ZonotopeEnclosure}
  max_steps: Int64 10000
  abs_tol: Float64 1.0e-10
  orderT: Int64 8
  orderQ: Int64 2
  disjointness: ReachabilityAnalysis.ZonotopeEnclosure ReachabilityAnalysis.ZonotopeEnclosure()
  adaptive: Bool true
  min_abs_tol: Float64 1.0e-29

[mforets]

This problem turns out to be quite demanding for the solver, as you can see by taking a smaller time span and a smaller set of initial conditions:

julia> prob = @ivp(x' = cartpole!(x), dim: 4, x(0) ∈ X₀*0.01); # shrink the initial conditions

julia> @time sol = solve(prob, T=1.0, alg=TMJets(max_steps=10_000, abs_tol=1e-10));
 23.713397 seconds (229.12 M allocations: 15.720 GiB, 16.90% gc time)

julia> tspan(sol)
Interval(0.0, 1.0)

The full time horizon is met, though the runtime and allocations are quite high.

@SebastianGuadalupe can you see how much do we gain by rewriting cartpole! to better use @taylorize? (see here)

[lbenet]

See also this for some tips (and limitations) when using @taylorize; my guess is that you are not exploiting all its capacity. Admittedly, the code becomes not so nice...

[mforets]

With the auxiliary variable aux = (f + m*l*u[4]^2*sin(u[3])) (appears 3 times) the runtime reduces by more than 2x:

@time sol = solve(prob, T=1.0, alg=TMJets(max_steps=10_000, abs_tol=1e-10));
  8.590104 seconds (175.30 M allocations: 12.130 GiB, 23.26% gc time)

[mforets]

PR #214 let's one pass these options from outside the struct TMJets. (if both the stuct and the options are passed, those in the struct prevail). i'll add some comments in the documentation about this behavior later.