decomposed ASB07 algorithm

src/ReachSets/ContinuousPost/ASB07_decomposed/ASB07_decomposed.jl

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

src/ReachSets/ContinuousPost/ASB07_decomposed/init.jl

@@ -0,0 +1,16 @@
+init(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒜, 𝑆, copy(𝑂))
+
+function init!(𝒜::ASB07_decomposed, 𝑆::AbstractSystem, 𝑂::Options)
+    𝑂copy = copy(𝑂)
+    𝑂copy[:n] = statedim(𝑆)
+
+    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)
+
+    𝑂validated[:partition] = 𝒜.options[:partition]
+
+    # compute blocks
+    𝑂validated[:blocks] =
+        compute_blocks(𝒜.options[:vars], 𝑂validated[:partition])
+
+    return 𝑂validated
+end

src/ReachSets/ContinuousPost/ASB07_decomposed/options.jl

@@ -0,0 +1,30 @@
+function options_ASB07_decomposed()
+    𝑂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"))
+
+    push!(𝑂spec, OptionSpec(:block_options_init, Hyperrectangle, domain=Any,
+        info="option for the decomposition of the initial states"))
+
+    push!(𝑂spec, OptionSpec(:partition, [Int[]],
+        domain=AbstractVector{<:AbstractVector{Int}}, domain_check=ispartition,
+        info="block partition; a block is represented by a vector containing " *
+             "its indices"))
+
+    push!(𝑂spec, OptionSpec(:vars, Int[], domain=AbstractVector{Int},
+        domain_check=(v  ->  length(v) > 0 && all(e -> e > 0, v)),
+        info="variables of interest; default: all variables"))
+
+    return 𝑂spec
+end

src/ReachSets/ContinuousPost/ASB07_decomposed/post.jl

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

src/ReachSets/ContinuousPost/ASB07_decomposed/reach.jl

@@ -0,0 +1,58 @@
+function reach_ASB07_decomposed!(R::Vector{<:ReachSet},
+                                 Ωhat0::Vector{<:LazySet},
+                                 U::Union{ConstantInput, Nothing},
+                                 ϕ::AbstractMatrix{IntervalArithmetic.Interval{NUM}},
+                                 N::Int,
+                                 δ::Float64,
+                                 max_order::Int,
+                                 n::Int,
+                                 partition,
+                                 blocks) where {NUM}
+    # initial reach set
+    t0, t1 = zero(δ), δ
+    R[1] = ReachSet(CartesianProductArray(Ωhat0[blocks]), t0, t1)
+
+    b = length(blocks)
+    Rₖ_array = Vector{Zonotope{NUM}}(undef, b)
+    ϕpowerk = copy(ϕ)
+
+    if U != nothing
+        Whatk = Vector{Zonotope{NUM}}(undef, b)
+        inputs = next_set(U)
+        @inbounds for i in 1:b
+            bi = partition[blocks[i]]
+            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
+        end
+    end
+
+    k = 2
+    while k <= N
+        for i in 1:b
+            bi = partition[blocks[i]]
+            Rₖ_bi = ZeroSet(length(bi))
+            for (j, bj) in enumerate(partition)
+                block = IntervalMatrix(ϕpowerk[bi, bj])
+                if !iszero(block)
+                    Rₖ_bij = overapproximate(block * Ωhat0[j], Zonotope)
+                    Rₖ_bi = minkowski_sum(Rₖ_bi, Rₖ_bij)
+                end
+            end
+            if U != nothing
+                Rₖ_bi = minkowski_sum(Rₖ_bi, Whatk[i])
+            end
+            Rₖ_bi = reduce_order(Rₖ_bi, max_order)  # reduce order
+            Rₖ_array[i] = Rₖ_bi
+        end
+        Rₖ = CartesianProductArray(copy(Rₖ_array))
+
+        # store reach set
+        t0 = t1
+        t1 += δ
+        R[k] = ReachSet(Rₖ, t0, t1)
+
+        ϕpowerk *= ϕ
+
+        k += 1
+    end
+    return R
+end

src/ReachSets/ReachSets.jl

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