A brief introduction to JuliaReach

Marcelo Forets, Universidad de la República, Uruguay.

Berlin Julia Users Meetup 2019, Berlin, August 7th, 2019

Outline of this talk

  • Context of the JuliaReach org

  • LazySets: from low dimensions to high dimensions using multiple dispatch

  • Applications

    • Neural network verification
    • Reachability analysis of uncertain dynamical systems

What is JuliaReach?

  • A toolbox to do Reachability Computations for Dynamical Systems in Julia

What are dynamical systems?

Just about anything that evolves with time!

  • In particular we focus on cyber-physical systems: a device that by means of computation is able to control or interact with a physical process

  • Examples:

    • embedded controllers (aircrafts, autonomous cars, ...)
    • medical devices; safety-critical or mission-critical devices; power systems stability; ...

Who makes JuliaReach?

Bogomolov, S., Forets, M., Frehse, G., Potomkin, K., & Schilling, C. (2019, April). JuliaReach: a toolbox for set-based reachability. In Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control (pp. 39-44), ACM.

An intro to LazySets

In [10]:
using LazySets, Plots
using AbstractTrees

AbstractTrees.children(x::Type) = subtypes(x)
print_tree(LazySet)
LazySet
├─ AbstractCentrallySymmetric
│  ├─ Ball2
│  ├─ Ballp
│  └─ Ellipsoid
├─ AbstractPolyhedron
│  ├─ AbstractPolytope
│  │  ├─ AbstractCentrallySymmetricPolytope
│  │  │  ├─ AbstractZonotope
│  │  │  │  ├─ AbstractHyperrectangle
│  │  │  │  │  ├─ AbstractSingleton
│  │  │  │  │  │  ├─ Singleton
│  │  │  │  │  │  └─ ZeroSet
│  │  │  │  │  ├─ BallInf
│  │  │  │  │  ├─ Hyperrectangle
│  │  │  │  │  ├─ Interval
│  │  │  │  │  └─ SymmetricIntervalHull
│  │  │  │  ├─ LineSegment
│  │  │  │  └─ Zonotope
│  │  │  └─ Ball1
│  │  ├─ AbstractPolygon
│  │  │  ├─ AbstractHPolygon
│  │  │  │  ├─ HPolygon
│  │  │  │  └─ HPolygonOpt
│  │  │  └─ VPolygon
│  │  ├─ HPolytope
│  │  └─ VPolytope
│  ├─ HPolyhedron
│  ├─ HalfSpace
│  ├─ Hyperplane
│  ├─ Line
│  └─ Universe
├─ AffineMap
├─ CacheMinkowskiSum
├─ CartesianProduct
├─ CartesianProductArray
├─ ConvexHull
├─ ConvexHullArray
├─ EmptySet
├─ ExponentialMap
├─ ExponentialProjectionMap
├─ Intersection
├─ IntersectionArray
├─ LinearMap
├─ MinkowskiSum
├─ MinkowskiSumArray
├─ ResetMap
└─ Translation

Convex sets

Balls

More balls

Unbounded sets

Translation

Minkowski sum

Linear map

Projection

Union (not convex)

Convex hull

Intersection

A high-dimensional demo

  • We'l work with a 100-dimensional set which has... exp(100 * log(2)) ≈ 1.2676506002282317e30 vertices.
In [70]:
H = rand(Hyperrectangle, dim=100);
d = ones(Float64, 100);
In [71]:
@btime σ($d, $H);
  97.050 ns (1 allocation: 896 bytes)
In [72]:
σ(d, H)  H
Out[72]:
true
In [73]:
M = rand(100, 100)
P = project(M * H, [99, 100], LinearMap);
In [74]:
@btime project($M * $H, [99, 100], LinearMap);
  651.521 ns (15 allocations: 2.63 KiB)
In [75]:
plot(overapproximate(P, HPolygon, 1e-2))
Out[75]:
-30 -20 -10 0 10 20 -20 -10 0 10 20 30 40
In [80]:
@btime overapproximate($P, HPolygon, 1e-2)
  16.663 ms (10893 allocations: 2.97 MiB)
Out[80]:
HPolygon{Float64}(HalfSpace{Float64,VN} where VN<:AbstractArray{Float64,1}[HalfSpace{Float64,Array{Float64,1}}([1.0, 0.0], 59.3398), HalfSpace{Float64,Array{Float64,1}}([0.980785, 0.19509], 71.9015), HalfSpace{Float64,Array{Float64,1}}([0.92388, 0.382683], 81.7768), HalfSpace{Float64,Array{Float64,1}}([0.83147, 0.55557], 88.5864), HalfSpace{Float64,Array{Float64,1}}([0.707107, 0.707107], 92.0686), HalfSpace{Float64,Array{Float64,1}}([0.55557, 0.83147], 92.0894), HalfSpace{Float64,Array{Float64,1}}([0.382683, 0.92388], 88.6482), HalfSpace{Float64,Array{Float64,1}}([0.19509, 0.980785], 81.8771), HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0], 72.0363), HalfSpace{Float64,Array{Float64,1}}([-0.0444605, 0.999011], 69.4208)  …  HalfSpace{Float64,Array{Float64,1}}([0.964989, -0.262291], 39.9318), HalfSpace{Float64,Array{Float64,1}}([0.970962, -0.239233], 41.6829), HalfSpace{Float64,Array{Float64,1}}([0.977136, -0.212617], 43.7054), HalfSpace{Float64,Array{Float64,1}}([0.980848, -0.194775], 45.0502), HalfSpace{Float64,Array{Float64,1}}([0.983167, -0.182712], 45.9559), HalfSpace{Float64,Array{Float64,1}}([0.987516, -0.157518], 47.8345), HalfSpace{Float64,Array{Float64,1}}([0.98866, -0.150174], 48.3841), HalfSpace{Float64,Array{Float64,1}}([0.99458, -0.103978], 51.8667), HalfSpace{Float64,Array{Float64,1}}([0.999056, -0.0434322], 56.2681), HalfSpace{Float64,Array{Float64,1}}([0.999633, -0.0271011], 57.4238)])
In [81]:
f(B) = M * H  B
B = Ball2(100 * rand(100), 2.0)
P = project(f(B), [99, 100], LinearMap);
In [79]:
@btime f($B)
  61.150 ns (2 allocations: 64 bytes)
Out[79]:
MinkowskiSum{Float64,LinearMap{Float64,Hyperrectangle{Float64},Float64,Array{Float64,2}},Ball2{Float64}}(LinearMap{Float64,Hyperrectangle{Float64},Float64,Array{Float64,2}}([0.850733 0.829261 … 0.24797 0.262343; 0.589077 0.342921 … 0.381546 0.0782568; … ; 0.300484 0.279613 … 0.0927591 0.446477; 0.909262 0.294161 … 0.494402 0.734955], Hyperrectangle{Float64}([-0.703713, 0.157051, -0.672884, -0.254412, 1.07036, 0.157169, -0.0758222, 0.218267, -2.83684, -0.0165159  …  1.16904, 0.710912, 0.35657, 1.61172, -1.20456, -0.41658, -0.197501, 1.65477, 2.71736, 1.25608], [0.137156, 0.496576, 0.122476, 0.33142, 2.2767, 0.622647, 0.00417522, 0.0716649, 1.50485, 1.17094  …  0.193609, 1.34673, 1.05655, 0.138119, 1.31267, 0.0610448, 0.310624, 2.00374, 0.0733485, 1.04473])), Ball2{Float64}([76.4493, 92.8461, 91.3955, 7.13363, 17.7478, 20.4587, 29.317, 33.3211, 18.9079, 43.6403  …  49.3594, 23.784, 5.7038, 12.9211, 72.8995, 56.2226, 92.5188, 87.691, 29.8813, 31.076], 2.0))
In [68]:
plot!(overapproximate(P, HPolygon, 1e-3))
Out[68]:
-40 -20 0 20 40 -25 0 25 50 75 100
In [82]:
@btime overapproximate($P, HPolygon, 1e-3)
  37.953 ms (24945 allocations: 6.80 MiB)
Out[82]:
HPolygon{Float64}(HalfSpace{Float64,VN} where VN<:AbstractArray{Float64,1}[HalfSpace{Float64,Array{Float64,1}}([1.0, 0.0], 67.2004), HalfSpace{Float64,Array{Float64,1}}([0.998795, 0.0490677], 71.4796), HalfSpace{Float64,Array{Float64,1}}([0.995185, 0.0980171], 75.5914), HalfSpace{Float64,Array{Float64,1}}([0.989177, 0.14673], 79.526), HalfSpace{Float64,Array{Float64,1}}([0.980785, 0.19509], 83.2737), HalfSpace{Float64,Array{Float64,1}}([0.970031, 0.24298], 86.8257), HalfSpace{Float64,Array{Float64,1}}([0.95694, 0.290285], 90.1734), HalfSpace{Float64,Array{Float64,1}}([0.941544, 0.33689], 93.3086), HalfSpace{Float64,Array{Float64,1}}([0.92388, 0.382683], 96.2238), HalfSpace{Float64,Array{Float64,1}}([0.903989, 0.427555], 98.9121)  …  HalfSpace{Float64,Array{Float64,1}}([0.98866, -0.150174], 53.336), HalfSpace{Float64,Array{Float64,1}}([0.98972, -0.143019], 54.0246), HalfSpace{Float64,Array{Float64,1}}([0.99458, -0.103978], 57.7324), HalfSpace{Float64,Array{Float64,1}}([0.995226, -0.0975957], 58.329), HalfSpace{Float64,Array{Float64,1}}([0.99592, -0.0902353], 59.0165), HalfSpace{Float64,Array{Float64,1}}([0.997866, -0.0652955], 61.3197), HalfSpace{Float64,Array{Float64,1}}([0.999056, -0.0434322], 63.3058), HalfSpace{Float64,Array{Float64,1}}([0.999548, -0.0300498], 64.5061), HalfSpace{Float64,Array{Float64,1}}([0.999633, -0.0271011], 64.7726), HalfSpace{Float64,Array{Float64,1}}([0.999896, -0.0144429], 65.9109)])

Fair enough, show me the trick.

[KAR19] Stefan Karpinksi (Julia Computing). The Unreasonable Effectiveness of Multiple Dispatch. JuliaCon 2019

https://www.youtube.com/watch?v=kc9HwsxE1OY

The lazy evaluation paradigm

  • Principle of lazy evaluation: delay the evaluation of an expression until its value is needed.
  • Question: what do i really need to compute if i just want $\rho(d, MH \oplus B)$?
  • Answer: combine basic machinery in convex geometry with the lazy evaluation paradigm, letting Julia to handle the appropriate dispatch for your computation.

Calculus with suport functions

Calculus with suport functions (cont.)

Case study: lazy Minkowski sum

  • Implement the Minkowski sum operator as a new LazySet, with two basic functions:
    • dimension
    • compute its support vector (or support function)
function ρ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}
    return ρ(d, ms.X) + ρ(d, ms.Y)
end

Hold on, did he just say that is a set?

struct MinkowskiSum{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2

    # default constructor with dimension check
    function MinkowskiSum{N, S1, S2}(X::S1, Y::S2) where
            {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a Minkowski sum must have the " *
            "same dimension"
        return new{N, S1, S2}(X, Y)
    end
end
@neutral(MinkowskiSum, ZeroSet)
@absorbing(MinkowskiSum, EmptySet)
(X::LazySet, Y::LazySet) = MinkowskiSum(X, Y)

Reachability Analysis in JuliaReach

The problem

  • Complex, real-world systems are prone to failures
  • New industrial standars explicitly recommend the user of formal methods: DO-178C (avionics), ISO 26262 (automotive systems), IEC 62304 (medical devices), and CENELEC EN 50128 (railway systems).

  • Engineers need safety and reliability guarantees to take design decisions under non-deterministic inputs, parameters or noise

$\Longrightarrow$ The field of reachability analysis is concerned with understanding the set of all possible behaviors of such systems.

Bogomolov, S., Forets, M., Frehse, G., Potomkin, K., & Schilling, C. (2019, April). JuliaReach: a toolbox for set-based reachability. In Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control (pp. 39-44), ACM.

References to our research work

  • [BFFPSV18] Reach Set Approximation through Decomposition with Low-dimensional Sets and High-dimensional Matrices. S. Bogomolov, M. Forets, G. Frehse, , A. Podelski, C. Schilling and F. Viry (2018) HSCC'18 Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control: 41–50.

  • [ARCH-COMP18-AFF] ARCH-COMP18 Category Report: Continuous and Hybrid Systems with Linear Continuous Dynamics. M. Althoff, S. Bak, X. Chen, C. Fan, M. Forets, G. Frehse, N. Kochdumper, Y. Li, S. Mitra, R. Ray, C. Schilling and S. Schupp (2018) ARCH18. 5th International Workshop on Applied Verification of Continuous and Hybrid Systems, 54: 23–52. doi: 10.29007/73mb.

  • [BFFPS19] JuliaReach: a Toolbox for Set-Based Reachability. Sergiy Bogomolov, Marcelo Forets, Goran Frehse, Kostiantyn Potomkin, Christian Schilling. Published in Proceedings of HSCC'19: 22nd ACM International Conference on Hybrid Systems: Computation and Control (HSCC'19).

References to our research work (cont.)

  • [BFFPS19b] Reachability analysis of linear hybrid systems via block decomposition. Sergiy Bogomolov, Marcelo Forets, Goran Frehse, Kostiantyn Potomkin, Christian Schilling.

  • [ARCH-COMP19-AFF] ARCH-COMP19 Category Report: Continuous and Hybrid Systems with Linear Continuous Dynamics. Matthias Althoff, Stanley Bak, Marcelo Forets, Goran Frehse, Niklas Kochdumper, Rajarshi Ray, Christian Schilling and Stefan Schupp (2019) ARCH19. 6th International Workshop on Applied Verification of Continuous and Hybrid Systems, vol 61, pages 14–40 doi: 10.29007/bj1w.

  • [ARCH-COMP19-NLN] ARCH-COMP19 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics. Fabian Immler, Matthias Althoff, Luis Benet, Alexandre Chapoutot, Xin Chen, Marcelo Forets, Luca Geretti, Niklas Kochdumper, David P. Sanders and Christian Schilling (2019) ARCH19. 6th International Workshop on Applied Verification of Continuous and Hybrid Systems, vol 61, pages 41–61 doi: 10.29007/bj1w.

ISS Benchmark

  • ISS: structural model of component 1R (Russian service module) of the international space station
  • 270 state variables with three uncertain inputs
  • linear time-invariant: $x'(t) = Ax(t) + Bu(t)$, $u(t) \in \mathcal{U}$, $y(t) = Cx(t) + Du(t)$
  • specification: output $y_3(t)$ is within the range $[0.0007, 0.0007]$ for all $t \in [0, 20]$.

ISS Benchmark

ISS Benchmark

ISS Benchmark

Quadrotor benchmark

  • nonlinear
  • 12 dimensions
  • safety property: stabilization of the height after some given time horizon, which should stay within given bounds
  • uncertain positions and uncertain velocities

Quadrotor benchmark

Quadrotor benchmark

Quadrotor benchmark

Quadrotor benchmark

Verifying Deep Neural Networks

Any sufficiently advanced technology is indistinguishable from magic.

(Clark's third law)

Conclusions

  • Convex sets are expressive, calculus is tractable.
  • Closure under most standard set operations.
  • Support functions allow for lazy (on-demand) computations.
  • Non-convex sets: approximate by (union of) convex sets.
  • Julia is an excellent choice to reason about set-based algorithms.

Acknowledgements

  • JuliaReach development is led by:

    • Marcelo Forets (Univ. de la República, Uruguay)
    • Christian Schilling (Institute of Science and Technology, Austria)
  • Contributors to JuliaReach open projects include:

    • A. Deshmuhk (Indian Institute of Information Technology, India)
    • B. Garate (Univ. de la República, Uruguay)
    • S. Guadalupe (Univ. de la República, Uruguay)
    • N. Kekatos (Univ. Grenoble Alpes, France)
    • B. Legat (UCLouvain, Belgium)
    • K. Potomkin (ANU, Australia)
    • A. Rocca (INRIA, France)
    • F. Viry (CERFACS, France)
  • Scientific collaborators:

    • L. Benet (UNAM, México)
    • S. Bogomolov (ANU, Australia)
    • G. Frehse (ENSTA ParisTech, France)
    • A. Podelski (Univ. of Freiburg, Germany),
    • David P. Sanders (UNAM, México) </font>

Join the effort! https://github.com/juliareach
<https://gitter.im/JuliaReach/Lobby

Contact: [email protected] http://github.com/mforets