Controller dependent reachable set for continuous ODEs with non-linear dynamics

[yuzhou42]

Hi, thank you for your work! May I know if it is possible to compute a controller-dependent reachable set for continuous ODEs with non-linear dynamics? I checked out the example for Quadrotor, but it only supports one fixed target. How can I get the reachable set with a set of control inputs or a certain control policy? Thank you!

[mforets]

Hi @yuzhou42, thanks for your interest in our work. We don't have an out-of-the-box user/solver interface for the type of problem that you want to solve: some manual work has to be done to prepare and solve such problems using JuliaReach.

Can you give us more details on the problem you want to solve, or a solution that you have attempted? Or an article doing a similar work? We would be happy to help if you have problems using our library.

Let me also comment on the other repo, https://github.com/JuliaReach/NeuralNetworkAnalysis.jl . In that case, the controls are given by neural network outputs, but if you have say, a Julia function instead, i think it shouldn't be hard to adapt such solver to use it. In fact you could check the simple black_box_toy_model.jl; the controller in this case is the function called "oscillator" and is state-dependent and it is applied at each multiple of the period.

[yuzhou42]

Hi @mforets,

Thank you for your reply! I'm facing two issues. First is that I'm not sure how to feed in the set of admissible control input signals for a continuous ODEs as there are no similar examples. Second is that I tried the cartpole model with ReachabilityAnalysis, but couldn't get it run for more than 1s. Here is the sample code:

@taylorize function cartpole!(dx, x, params, t)
    local m_c, m_p, l, g, u = 0.85944, 17.08072, 0.4626, 9.81, 0
    dx[1] = x[3]
    dx[2] = x[4]
    dx[3] = 1.0/(m_c+m_p*sin(x[2])^2)*(u+m_p*sin(x[2])*(l*x[4]^2-g*cos(x[2])))
    dx[4] = 1.0/(l*(m_c+m_p*sin(x[2])^2))*(-u*cos(x[2])-m_p*l*x[4]^2*cos(x[2])*sin(x[2])+(m_c + m_p)*g*sin(x[2]))
    return dx
end
X0 = Hyperrectangle(low=[0,-1,0,0], high=[0,1,0,0]);
prob = @ivp(x' = cartpole!(x), dim=4, x(0) ∈ X0 * 0.001)
@time sol = solve(prob, tspan=(0.0,0.5), alg=TMJets(maxsteps=1000000, abstol=1e-10));

The error output is like

┌ Warning: Minimum absolute tolerance reached: │ t[0] = 0.8112222762549506 │ E′ = IntervalArithmetic.Interval{Float64}[Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf)] │ E = IntervalArithmetic.Interval{Float64}[Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf)] │ _success = false │ all(iscontractive.(E′, E)) = false │ reduced_abstol = 1.0000000000000003e-51 └ @ TaylorModels /home/yu/.julia/packages/TaylorModels/b295K/src/validatedODEs.jl:747 ┌ Warning: Maximum number of validate steps reached. │ t[0] = 0.8112222762549506 │ E′ = IntervalArithmetic.Interval{Float64}[Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf)] │ E = IntervalArithmetic.Interval{Float64}[Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf), Interval(-Inf, Inf)] │ _success = false │ all(iscontractive.(E′, E)) = false └ @ TaylorModels /home/yu/.julia/packages/TaylorModels/b295K/src/validatedODEs.jl:762 ┌ Warning: Exiting due to unsuccessfull step │ _success = false │ t0 = 0.8112236622119696 └ @ TaylorModels /home/yu/.julia/packages/TaylorModels/b295K/src/validatedODEs.jl:894

I checked out the black_box_toy_model.jl example you pointed out, it is somewhat similar to what I wanted to do. But unfortunately, I couldn't add NeuralNetworkAnalysis successfully to run it.

Best, Yu

[mforets]

. But unfortunately, I couldn't add NeuralNetworkAnalysis successfully to run it.

just in case: the package is not registered yet so you need to do

] add https://github.com/JuliaReach/NeuralNetworkAnalysis.jl.git

or did you get an error elsewhere?

Concerning the cartpole model until 1s, it may help to start with subsets of the initial states and find algorithm parameters that give reasonable results (in terms of runtime and accuracy).

Moreover, it is possible to get a factor of ~2x perf improvement if you follow the tips about the @taylorize macro explained here, because that helps the way the macro works. In this case I wrote:

@taylorize function cartpole_v2!(dx, x, params, t)
    local m_c, m_p, l, g, u = 0.85944, 17.08072, 0.4626, 9.81, 0
    
    aux0 = sin(x[2])
    aux1 = aux0^2
    aux3 = cos(x[2])
    aux4 = x[4]^2
    aux5 = m_c + m_p*aux1
    aux6 = l*aux4 - g*aux3
    aux7 = m_p*l*aux4*aux3*aux0
    aux8 = (m_c + m_p)*g*aux0
    aux9 = u + m_p*aux0*aux6
    aux10 = -u*aux3-aux7+aux8
    aux11 = 1.0 / aux5
    aux12 = aux10 / l

    dx[1] = x[3]
    dx[2] = x[4]
    dx[3] = aux11 * aux9
    dx[4] = aux11 * aux12
    return dx
end

Here are some results. I didn't play too much with parameters, so I expect there is room for improvement. Moreover, it may be worth to try the multi-threaded flowpipe computation (start julia with a number of threads > 1, the default threading option should dispatch correctly): I'd be interested to know what results do you get.

X0 = Hyperrectangle(low=[0, -0.001,0,0], high=[0, 0.001,0,0]);

# here i chose a deliberately large number of splits in the x2 direction.. try with smaller values!
X0s = split(X0, [1, 100, 1, 1])

prob = @ivp(x' = cartpole_v2!(x), dim=4, x(0) ∈ X0s);

alg = TMJets21a(maxsteps=10_000, abstol=1e-5, orderQ=1, orderT=4, adaptive=true)
@time sol = solve(prob, tspan=(0.0, 1.0), alg=alg);
345.060740 seconds (3.74 G allocations: 325.414 GiB, 18.44% gc time)

# plots (long time)
#plot(sol, vars=(0, 1), xlab="t", ylab="x₁", alpha=0.8, lw=0.0, c=:blue)
#plot(sol, vars=(0, 2), xlab="t", ylab="x₂", alpha=0.8, lw=0.0, c=:blue)

[yuzhou42]

Hi @mforets, thank you for your help. With multithreading, the calculation time for the same setup is about 126s on my side. Decreasing the number of splits to 10 also speeds up the time to 12s for 1s reachable set. But I find it's still a bit hard to tune when I increase the time for the reachable set to more than 1s. May I know how I can add constraints on the states? Thanks!

[mforets]

With multithreading, the calculation time for the same setup is about 126s on my side. Decreasing the number of splits to 10 also speeds up the time to 12s for 1s reachable set.

:+1: Good to know that with 10 splits you can reach to the desired time horizon. As usual, there is a tradeoff between the quality of the approximation and the computational resources used, and with more splitting to cover the initial set you may expect smaller sets at intermediate and final time (there are exceptions though).

But I find it's still a bit hard to tune when I increase the time for the reachable set to more than 1s.

Unfortunately, the current state of affairs is that parameter tuning cannot be avoided in most scenarios. In fact, finding good parameter heuristics at runtime is a research problem, see e.g. [1].

May I know how I can add constraints on the states?

IIUC you'd like to evaluate one or more reach-sets, make some transformation on the state variables and continue the computation with those new values. An example of a computation along those lines can be found here, where X is the current reach-set in n + m variables (n state variable and m control variables), and U0 is a set that corresponds to the new control values; it is a set in R^m.

Understanding this would require familiarity with TaylorModels.jl, TaylorSeries.jl and IntervalArithmetic.jl which make the basis for the internal representation of Taylor model reach-sets used in this library, and also the range approximation strategies found in LazySets.jl and in ReachabilityAnalysis.jl of course.

We plan to add new publicly available documentation soon, and we are preparing a workshop to give at the JuliaCon 2021 conference on these aspects [2]. You are invited to attend the workshop!

PS feel free to reach us on gitter or the (JuliaLang) zulip channel reachability-analysis

[1] https://ieeexplore.ieee.org/document/9304431

[2] https://github.com/JuliaReach/JuliaCon-2021-Workshop-Its-All-Set

[mforets]

May I know how I can add constraints on the states?

Maybe I misunderstood your question. Do you have an example of the type of constraints that you want to add? For instance if you want to add conditions (spatial or temporal), that is achieved by adding an invariant to the system definition, like @ivp(x' = cartpole_v2!(x), dim=4, x(0) ∈ X0s, x ∈ W), where W is an n-dimensional set. Then, whenever a reach-set is disjoint with the invariant W, the simulation halts.

[yuzhou42]

Thank you for your help @mforets! You've solved all the questions I have so far. Looking forward to the new documentation and the workshop!