Week #4: Reachability Analysis

Summary of contributions:

• LazySets#2197

• NeuralVerification#11

• ReachabilityAnalysis#281

• ReachabilityAnalysis#283

• ReachabilityAnalysis#285

• ReachabilityAnalysis#286

• ReachabilityAnalysis#288

• ReachabilityAnalysis#289

• ReachabilityAnalysis#290

Introduction to Reachability Analysis

Reachability analysis is the study of which sets of states a dynamical system can reach, starting from a set of initial states and under the influence of a set of input trajectories and parameter values. These techniques have been successfully applied in the literature in diverse domains such as:

• Robustness verification: Determining whether a system is safe, i.e. to verify if it does not enter into a region of unsafe sets. Typical applications are assessing the critical distance between autonomous vehicles or robots, or critical concentration of chemicals in a reactor.

• Validation of control strategies: Checking if the system trajectories stay in a region around a reference trajectory, or reach a goal region around a setpoint, for any admissible value of the non-deterministic inputs, initial conditions or noise.

• Controller synthesis: Finding parameter sets of controllers that satisfy safety or performance constraints.

The formal definition of the reachability problem is given below. Given an ordinary differential equation (ODE)

where inline_formula not implemented denotes the set of initial states, inline_formula not implemented denotes the input set, and inline_formula not implemented denotes the set of parameter values, the reachable set at time inline_formula not implemented is defined by

The reachable set associated to a continuous time interval inline_formula not implemented, also known as the flowpipe, is defined as:

This definition is naturally extended for the flowpipe over a finite time span inline_formula not implemented. In practice, the set inline_formula not implemented, and the flowpipe, cannot be obtained exactly. Reachability methods aim at computing suitable over-approximations (or under-approximations) of the flowpipe.

Among the different reachability methods, set-based reachability uses set representations such as zonotopes or other (e.g. support functions, ellipsoids, polytopes, taylor models) to overapproximate the flowpipe of a given system. In particular, zonotopes have been used extensively because linear maps and minkowski sums can be computed very efficiently [9].

The Julia library ReachabilityAnalysis.jl and the packages in the JuliaReach ecosystem have implemented reachability methods for linear, nonlinear and hybrid ODEs, this is the library that I will be working with and extending

Examples

Van der Pol

1 - Using ReachabilityAnalysis.jl I will be demonstrating how it is used with a well known differential equation, the van der Pol oscillator. This model is a non-conservative oscillator with non-linear damping. It evolves in time according to the second-order differential equation [1]:

where x is the position coordinate—which is a function of the time t, and μ is a scalar parameter indicating the nonlinearity and the strength of the damping.

If we take inline_formula not implemented then we get:

this is a form that can be easily translated to the ReachabilityAnalysis system

using ReachabilityAnalysis

@taylorize function vanderpol!(dx, x, params, t)

    local μ = 1.0

    dx[1] = x[2]

    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]

    return dx

end

Here we define the function vanderpol! that will be containing the dynamics of our system, the arguments of the function are always the same, then the local variable μ is defined for ease of change and readability. After that the differential equation is defined, with dx[1] being inline_formula not implemented and dx[2] being inline_formula not implemented, then the equations are typed normally with the variables inline_formula not implemented being x[i], or in this case inline_formula not implemented being x[1], x[2] respectively

We now can define the initial conditions, let's say that our initial conditions are the following:

X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])

prob = @ivp(x' = vanderpol!(x), dim=2, x(0) ∈ X0)

sol = solve(prob, T=7.0, alg=TMJets());

Here we defined the initial conditions X0(line 1), created the "Initial Value Problem" that ReachabilityAnalysis requires to solve the problem (line 2), and solved the problem for a time span of 7 seconds with the algorithm TMJets, that solves nonlinear models (line 3). TMJets solves problems using Taylor Models, for example Taylor polynomials with guaranteed error bounds to approximate functions.

For further computations, it is convenient to work with a zonotopic overapproximation of the flowpipe.

solz = overapproximate(sol, Zonotope);

using Plots

plot(solz, vars=(1, 2), lw=0.2, xlims=(-2.5, 2.5), xlab="x", ylab="y", label="Flowpipe", legend=:bottomright)

plot!(X0, label="X(0)")

In this plot we can see the initial condition and the resulting flowpipe, after this point, we could check if the flowpipe is compliant with a certain property.

Lotka-Volterra

2 - Lotka-Volterra equations: These equations also known as the predator–prey equations, are a pair of first-order nonlinear differential equations, frequently used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey. The populations change through time according to the pair of equations:

where inline_formula not implemented is the number of prey (for example, rabbits);

inline_formula not implemented is the number of some predator (for example, foxes);

inline_formula not implemented and inline_formula not implemented represent the instantaneous growth rates of the two populations;

inline_formula not implemented represents time;

inline_formula not implemented are positive real parameters describing the interaction of the two species.

@taylorize function lotkavolterra!(dx, x, params, t)

    local α, β, γ, δ = 2/3, 4/3, 1., 1.

    dx[1] = α * x[1] - β * x[1] * x[2]

    dx[2] = δ * x[1] * x[2] - γ * x[2]

    return dx

end

X0 = Hyperrectangle(low=[1., 1.], high=[1.1, 1.1])

prob = @ivp(x' = lotkavolterra!(x), dim=2, x(0) ∈ X0)

sol = solve(prob, T=8.0, alg=TMJets());

solz = overapproximate(sol, Zonotope);

using Plots

plot(solz, vars=(1, 2), alpha=0.3,lw=0., xlab="x", ylab="y", label="Flowpipe", legend=:bottomright)

plot!(X0, label="X(0)")

In this example I solved the differential equation for an initial condition inline_formula not implemented for a time span of 8 seconds using TMJets

Adding parameter variation

ReachabilityAnalysis allows you to add uncertainty in parameters, you can introduce them by adding a dimension to the problem and assigning it a derivative equal to zero

@taylorize function f(du, u, p, t)

    du[1] = u[3] * u[1] - u[4] * (u[1] * u[2]) - u[7] * u[1]^2

    du[2] = -u[5] * u[2] + u[6] * (u[1] * u[2])

    #encode uncertain params

    du[3] = zero(u[1]) # p[1]

    du[4] = zero(u[1]) # p[2]

    du[5] = zero(u[1]) # p[3]

    du[6] = zero(u[1]) # p[4]

    du[7] = zero(u[1]) # p[5]

end

# encode initial-value problem

p_int = (0.99..1.01) × (0.99..1.01) × (2.99..3.01) × (0.99..1.01) × (0.099..0.101)

U0 = Singleton([1.0, 1.0]) × convert(Hyperrectangle, p_int)

prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0);

sol = solve(prob, tspan=(0.0, 10.0));

solz = overapproximate(sol, Zonotope);

plot(solz, vars=(1, 2), lw=0.3, title="Uncertain params", lab="abs_tol = 1e-15", xlab="u1", ylab="u2")

Uncertain initial condition (inline_formula not implemented)

Now we consider an initial box around a point inline_formula not implemented with a radius of inline_formula not implemented

inline_formula not implemented

u0 = Singleton([1.0, 1.0])

□(ϵ) = BallInf(zeros(2), ϵ)

U0 = (u0 ⊕ □(0.05)) × convert(Hyperrectangle, p_int)

prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0)

sol = solve(prob, tspan=(0.0, 10.0), TMJets(abs_tol=1e-10))

solz = overapproximate(sol, Zonotope)

plot(solz, vars=(1, 2), color=:orange, lw=0.3,

     lab="eps = 0.05", title="Uncertain u0 and uncertain params",

     xlab="u1", ylab="u2")

inline_formula not implemented

U0 = (u0 ⊕ □(0.01)) × convert(Hyperrectangle, p_int)

prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0)

sol = solve(prob, tspan=(0.0, 10.0), TMJets(abs_tol=1e-10))

solz = overapproximate(sol, Zonotope)

plot!(solz, vars=(1, 2), color=:blue, lw=0.3,

  lab="eps = 0.01", title="Uncertain u0 and uncertain params",

  xlab="u1", ylab="u2")

Those were a couple of examples of the models that I have added to the documentation of ReachabilityAnalysis, now I will do a summary of the rest:

Other examples

Spacecraft: Spacecraft rendezvous is a perfect use case for formal verification of hybrid systems with nonlinear dynamics since mission failure can cost lives and is extremely expensive.

Platooning: The platooning benchmark considers a platoon of three vehicles following each other. This benchmark considers loss of communication between vehicles.

Building: A building benchmark from the ARCH competition.

SEIR model: The SEIR Model is an Compartmental model, these try to predict things such as how a disease spreads, or the total number infected, or the duration of an epidemic, and to estimate various epidemiological parameters such as the reproductive number.

Quadrotor: It studies the dynamics of a quadrotor.

Lorenz model: The model is a system of three ordinary differential equations now known as the Lorenz equations.

Future plans

In the following weeks I will be improving existing algorithms for hybrid ODEs using zonotopes, add a new reachability algorithm using zonotopes and I will expand the documentation of ReachabilityAnalysis.jl adding examples. Then I will be adding models to the package ReachabilityModels.jl, and setting up benchmarks in ReachabilityBenchamrks.jl for the testing of the different algorithms implemented in ReachabilityAnalysis.jl.

References

1. https://en.wikipedia.org/wiki/Van_der_Pol_oscillator

2. https://github.com/JuliaReach/ReachabilityAnalysis.jl

3. https://github.com/JuliaIntervals/TaylorModels.jl

4. https://cps-vo.org/group/ARCH/FriendlyCompetition