#639 - Reachability for LTI systems with interval matrices using zonotopes

src/ReachSets/ContinuousPost/ASB07/ASB07.jl

@@ -0,0 +1,21 @@
+export ASB07
+
+struct ASB07 <: ContinuousPost
+    options::TwoLayerOptions
+
+    function ASB07(𝑂::Options)
+        𝑂new = validate_and_wrap_options(𝑂, options_ASB07())
+        return new(𝑂new)
+    end
+end
+
+# convenience constructor from pairs of symbols
+ASB07(𝑂::Pair{Symbol,<:Any}...) = ASB07(Options(Dict{Symbol,Any}(𝑂)))
+
+# default options (they are added in the function validate_and_wrap_options)
+ASB07() = ASB07(Options())
+
+include("options.jl")
+include("init.jl")
+include("post.jl")
+include("reach.jl")

src/ReachSets/ContinuousPost/ASB07/init.jl

@@ -0,0 +1,7 @@
+init(𝒜::ASB07, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒜, 𝑆, copy(𝑂))
+
+function init!(::ASB07, 𝑆::AbstractSystem, 𝑂::Options)
+    𝑂[:n] = statedim(𝑆)
+
+    return validate_solver_options_and_add_default_values!(𝑂)
+end

src/ReachSets/ContinuousPost/ASB07/options.jl

@@ -0,0 +1,21 @@
+function options_ASB07()
+    𝑂spec = Vector{OptionSpec}()
+
+    push!(𝑂spec, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
+                            domain_check=(v  ->  v > 0.), info="time step"))
+
+    push!(𝑂spec, OptionSpec(:max_order, 10, domain=Int,
+                            info="maximum allowed order of zonotopes"))
+
+    push!(𝑂spec, OptionSpec(:order_discretization, 2, domain=Int,
+                            info="order of Taylor approximation in " *
+                                 "discretization"))
+
+    push!(𝑂spec, OptionSpec(:set_operations_discretization, "zonotope",
+                            domain=String,
+                            domain_check=(v  ->  v ∈ ["zonotope", "lazy"]),
+                            info="type of set operations applied during " *
+                            "discretization"))
+
+    return 𝑂spec
+end

src/ReachSets/ContinuousPost/ASB07/post.jl

@@ -0,0 +1,50 @@
+function post(𝒜::ASB07,
+              𝑃::InitialValueProblem{<:AbstractContinuousSystem},
+              𝑂::Options)
+    # =================================
+    # Initialization and discretization
+    # =================================
+
+    𝑂 = merge(𝒜.options.defaults, 𝑂, 𝒜.options.specified)
+    δ, T = 𝑂[:δ], 𝑂[:T]
+    N = round(Int, T / δ)
+
+    # compute and unrwap discretized system
+    𝑃_discrete = discretize(𝑃, δ; algorithm="interval_matrix",
+                            order=𝑂[:order_discretization],
+                            set_operations=𝑂[:set_operations_discretization])
+    Ω0, Φ = 𝑃_discrete.x0, 𝑃_discrete.s.A
+
+    # ====================
+    # Flowpipe computation
+    # ====================
+
+    # preallocate output
+    T = 𝑂[:set_operations_discretization] == "zonotope" ? Zonotope : LazySet
+    Rsets = Vector{ReachSet{T{Float64}}}(undef, N)
+
+    max_order = 𝑂[:max_order]
+
+    info("Reachable States Computation...")
+    @timing begin
+    if inputdim(𝑃_discrete) == 0
+        U = nothing
+    else
+        U = inputset(𝑃_discrete)
+    end
+    reach_ASB07!(Rsets, Ω0, U, Φ, N, δ, max_order)
+    end # timing
+
+    Rsol = ReachSolution(Rsets, 𝑂)
+
+    # ==========
+    # Projection
+    # ==========
+
+    if 𝑂[:project_reachset]
+        info("Projection...")
+        Rsol = @timing project(Rsol)
+    end
+
+    return Rsol
+end

src/ReachSets/ContinuousPost/ASB07/reach.jl

@@ -0,0 +1,28 @@
+function reach_ASB07!(R::Vector{<:ReachSet},
+                      Ω0::LazySet,
+                      U::Union{ConstantInput, Nothing},
+                      Φ::AbstractMatrix,
+                      N::Int,
+                      δ::Float64,
+                      max_order::Int)
+    # initial reach set
+    t0, t1 = zero(δ), δ
+    R[1] = ReachSet(Ω0, t0, t1)
+
+    k = 2
+    while k <= N
+        Rₖ = overapproximate(Φ * R[k-1].X, Zonotope)
+        if U != nothing
+            Rₖ = minkowski_sum(Rₖ, next_set(U))
+        end
+        Rₖ = reduce_order(Rₖ, max_order)  # reduce order
+
+        # store reach set
+        t0 = t1
+        t1 += δ
+        R[k] = ReachSet(Rₖ, t0, t1)
+
+        k += 1
+    end
+    return R
+end

src/ReachSets/ContinuousPost/GLGM06/init.jl

@@ -1,5 +1,5 @@
 # out-of-place initialization
-init(𝒫::GLGM06, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂))
+init(𝒜::GLGM06, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒜, 𝑆, copy(𝑂))
 
 function options_GLGM06()
 
@@ -8,11 +8,11 @@ function options_GLGM06()
     # step size
     push!(𝑂spec, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
                             domain_check=(v  ->  v > 0.), info="time step"))
- 
+
     # discretization
     push!(𝑂spec, OptionSpec(:discretization, "forward", domain=String,
                             info="model for bloating/continuous time analysis"))
-            
+
     push!(𝑂spec, OptionSpec(:sih_method, "concrete", domain=String,
                             info="method to compute the symmetric interval hull in discretization"))
 
@@ -27,7 +27,7 @@ function options_GLGM06()
 end
 
 # in-place initialization
-function init!(𝒫::GLGM06, 𝑆::AbstractSystem, 𝑂::Options)
+function init!(::GLGM06, 𝑆::AbstractSystem, 𝑂::Options)
 
     # state dimension
     𝑂[:n] = statedim(𝑆)

src/ReachSets/ContinuousPost/GLGM06/reach.jl

@@ -1,66 +1,64 @@
-# ===============================================================
+# ================
 # Homogeneous case
-# ===============================================================
-function reach_homog!(HR::Vector{ReachSet{Zonotope{Float64}}},
+# ================
+
+function reach_homog!(R::Vector{ReachSet{Zonotope{Float64}}},
                       Ω0::Zonotope,
                       Φ::AbstractMatrix,
                       N::Int,
                       δ::Float64,
                       max_order::Int)
-
-    # save timestamps with the reach set
-    t0, t1 = zero(δ), δ
-
     # initial reach set
-    HR[1] = ReachSet(Ω0, t0, t1)
+    t0, t1 = zero(δ), δ
+    R[1] = ReachSet(Ω0, t0, t1)
 
     k = 2
     while k <= N
-        HR_next = linear_map(Φ, HR[k-1].X)
-        HR_next_red = reduce_order(HR_next, max_order)
-        t0 = t1; t1 += δ
-        HR[k] = ReachSet(HR_next_red, t0, t1)
+        Rₖ = linear_map(Φ, R[k-1].X)
+        Rₖ = reduce_order(Rₖ, max_order)
+        t0 = t1
+        t1 += δ
+        R[k] = ReachSet(Rₖ, t0, t1)
+
         k += 1
     end
-    return HR
+    return R
 end
 
-# ===============================================================
+# ==================
 # Inhomogeneous case
-# ===============================================================
-function reach_inhomog!(X::Vector{ReachSet{Zonotope{Float64}}},
+# ==================
+
+function reach_inhomog!(R::Vector{ReachSet{Zonotope{Float64}}},
                         Ω0::Zonotope,
                         U::ConstantInput,
                         Φ::AbstractMatrix,
                         N::Int,
                         δ::Float64,
                         max_order::Int)
-
-    # initialization
+    # initial reach set
     t0, t1 = zero(δ), δ
-    V = next_set(U)
+    R[1] = ReachSet(Ω0, t0, t1)
 
-    X[1] = ReachSet(Ω0, t0, t1)
+    V = next_set(U)
     Wk₊ = V
     Φ_power_k = copy(Φ)
     Φ_power_k_cache = similar(Φ)
 
-    # update
     k = 2
     while k <= N
-        Xk = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊)
-        Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, V))
-
-        t0 = t1; t1 += δ
-        Xk_red = reduce_order(Xk, max_order)
-        X[k] = ReachSet(Xk_red, t0, t1)
+        Rₖ = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊)
+        Rₖ = reduce_order(Rₖ, max_order)
+        t0 = t1
+        t1 += δ
+        R[k] = ReachSet(Rₖ, t0, t1)
 
+        Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, V))
         Wk₊ = reduce_order(Wk₊, max_order)
 
         mul!(Φ_power_k_cache, Φ_power_k, Φ)
         copyto!(Φ_power_k, Φ_power_k_cache)
         k += 1
     end
-
-    return X 
+    return R
 end

src/ReachSets/ReachSets.jl

@@ -104,6 +104,8 @@ include("ContinuousPost/TMJets/TMJets.jl")
 
 include("ContinuousPost/BFFPS19/BFFPS19.jl")
 
+include("ContinuousPost/ASB07/ASB07.jl")
+
 # ========================
 # Reachability Algorithms
 # ========================

src/ReachSets/discretize.jl

@@ -801,7 +801,8 @@ function discretize_interval_matrix(𝑆::InitialValueProblem, δ::Float64,
     linear_maps = Vector{LinearMap{N}}(undef, order > 2 ? 3 : 2)
 
     A² = A * A
-    IδW = δ*I + 1/2 * δ^2 * A + 1/6 * δ^3 * A²
+    Iδ = IntervalMatrix(Diagonal(fill(IntervalMatrices.Interval(δ), size(A, 1))))
+    IδW = Iδ + 1/2 * δ^2 * A + 1/6 * δ^3 * A²
     linear_maps[1] = LinearMap(IδW, U0)
 
     E = _expm_remainder(A, δ, order; n=n)
@@ -870,7 +871,7 @@ end
 function _discretize_interval_matrix_inhomog(U, Ω0_homog, linear_maps,
                                              set_ops::Val{:zonotope})
     Ω0_inhomog = overapproximate(linear_maps[1], Zonotope)
-    @inbounds for i in 2:length(linear_map)
+    @inbounds for i in 2:length(linear_maps)
         Z = overapproximate(linear_maps[i], Zonotope)
         Ω0_inhomog = minkowski_sum(Z, Ω0_inhomog)
     end

test/Reachability/solve_continuous.jl

@@ -181,9 +181,9 @@ s = solve(IVP(LCS(A), X0), :T=>Inf, :mode=>"check", :property=>property)
 s = solve(IVP(CLCCS(A, B, X, U), X0),
           :T=>Inf, :mode=>"check", :property=>property)
 
-# =================================================================
-# Affine ODE with a nondeterministic input, x' = Ax + b, no inputs
-# =================================================================
+# ==================================
+# Affine ODE, x' = Ax + b, no inputs
+# ==================================
 A = [ 0.0509836  0.168159  0.95246   0.33644
       0.42377    0.67972   0.129232  0.126662
       0.518654   0.981313  0.489854  0.588326
@@ -220,3 +220,32 @@ solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2));
 property=(t,x)->x[2] <= 2.75
 𝑂 = Options(:T=>7.0, :mode=>"check", :property=>property)
 solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2))
+
+# =====
+# ASB07
+# =====
+
+# example 1 from
+# "Reachability analysis of linear systems with uncertain parameters and inputs"
+
+# initial set
+X0 = BallInf([1.0, 1.0], 0.1)
+
+# linear ODE: x' = Ax
+A = IntervalMatrix([-1.0 ± 0.05 -4.0 ± 0.05;
+                    4.0 ± 0.05 -1.0 ± 0.05])
+P_lin = IVP(LCS(A), X0)
+
+# affine ODE: x' = Ax + Bu
+B = IntervalMatrix(hcat([1.0 ± 0.0; 1.0 ± 0.0]))
+U = ConstantInput(Interval(-0.05, 0.05))
+P_aff = IVP(CLCCS(A, B, nothing, U), X0)
+
+for P in [P_lin, P_aff]
+    # default options
+    s = solve(P, Options(:T => 0.1), op=ASB07())
+
+    # use specific options
+    opC = ASB07(:δ => 0.04, :max_order => 10, :order_discretization => 4)
+    s = solve(P, Options(:T => 5.0), op=opC)
+end