From aaa96a555723a92d4b1b382d8799073cc1eec7fe Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 16 Mar 2019 17:11:38 +0100 Subject: [PATCH 1/3] allow unbounded-time setting --- .../ContinuousPost/BFFPSV18/check_blocks.jl | 36 +++++----- .../ContinuousPost/BFFPSV18/check_property.jl | 12 +++- .../ContinuousPost/BFFPSV18/reach.jl | 69 +++++++++++++++---- .../ContinuousPost/BFFPSV18/reach_blocks.jl | 47 +++++++------ src/ReachSets/ReachSets.jl | 2 +- 5 files changed, 112 insertions(+), 54 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_blocks.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_blocks.jl index 33c42d84..c4e85dd7 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_blocks.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_blocks.jl @@ -29,12 +29,14 @@ function check_blocks(ϕ::SparseMatrixCSC{NUM, Int}, U::Union{ConstantInput, Nothing}, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, eager_checking::Bool, - prop::Property + prop::Property, + progress_meter::Union{Progress, Nothing} )::Int where {NUM} + update!(progress_meter, 1) violation_index = 0 if !check(prop, CartesianProductArray(Xhat0[blocks])) if eager_checking @@ -59,9 +61,8 @@ function check_blocks(ϕ::SparseMatrixCSC{NUM, Int}, end k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] Xhatk_bi = ZeroSet(length(bi)) @@ -106,12 +107,14 @@ function check_blocks(ϕ::AbstractMatrix{NUM}, U::Union{ConstantInput, Nothing}, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, eager_checking::Bool, - prop::Property + prop::Property, + progress_meter::Union{Progress, Nothing} )::Int where {NUM} + update!(progress_meter, 1) violation_index = 0 if !check(prop, CartesianProductArray(Xhat0[blocks])) if eager_checking @@ -138,9 +141,8 @@ function check_blocks(ϕ::AbstractMatrix{NUM}, arr_length = (U == nothing) ? length(partition) : length(partition) + 1 k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] arr = Vector{LazySet{NUM}}(undef, arr_length) @@ -187,12 +189,14 @@ function check_blocks(ϕ::SparseMatrixExp{NUM}, U::Union{ConstantInput, Nothing}, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, eager_checking::Bool, - prop::Property + prop::Property, + progress_meter::Union{Progress, Nothing} )::Int where {NUM} + update!(progress_meter, 1) violation_index = 0 if !check(prop, CartesianProductArray(Xhat0[blocks])) if eager_checking @@ -217,9 +221,8 @@ function check_blocks(ϕ::SparseMatrixExp{NUM}, end k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] ϕpowerk_πbi = row(ϕpowerk, bi) @@ -262,12 +265,14 @@ function check_blocks(ϕ::SparseMatrixExp{NUM}, U::Union{ConstantInput, Nothing}, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, eager_checking::Bool, - prop::Property + prop::Property, + progress_meter::Union{Progress, Nothing} )::Int where {NUM} + update!(progress_meter, 1) violation_index = 0 if !check(prop, CartesianProductArray(Xhat0[blocks])) if eager_checking @@ -293,9 +298,8 @@ function check_blocks(ϕ::SparseMatrixExp{NUM}, arr_length = (U == nothing) ? length(partition) : length(partition) + 1 k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] arr = Vector{LazySet{NUM}}(undef, arr_length) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index 0fb9aaef..3cfe0a75 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -35,7 +35,8 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, n = statedim(S) blocks = options[:blocks] partition = convert_partition(options[:partition]) - N = ceil(Int, options[:T] / options[:δ]) + T = options[:T] + N = (T == Inf) ? nothing : ceil(Int, T / options[:δ]) # Cartesian decomposition of the initial set if length(partition) == 1 && length(partition[1]) == n && @@ -138,8 +139,15 @@ function check_property(S::IVP{<:AbstractDiscreteSystem}, # add property push!(args, property) - # call the adequate function with the given arguments list info("- Computing successors") + + # progress meter + progress_meter = (N != nothing) ? + Progress(N, 1, "Computing successors ") : + nothing + push!(args, progress_meter) + + # call the adequate function with the given arguments list answer = @timing available_algorithms_check[algorithm_backend]["func"](args...) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl index a115e599..1594dd84 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -4,6 +4,7 @@ using LazySets: CacheMinkowskiSum, isdisjoint import LazySets.Approximations: overapproximate +import ProgressMeter: update! """ reach(S, invariant, options) @@ -50,7 +51,8 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, n = statedim(S) blocks = options[:blocks] partition = convert_partition(options[:partition]) - N = ceil(Int, options[:T] / options[:δ]) + T = options[:T] + N = (T == Inf) ? nothing : ceil(Int, T / options[:δ]) # Cartesian decomposition of the initial set if length(partition) == 1 && length(partition[1]) == n && @@ -76,7 +78,7 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, else res_type = ReachSet{Hyperrectangle{NUM}, NUM} end - res = Vector{res_type}(undef, N) + res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N) # shortcut if only the initial set is required if N == 1 @@ -160,14 +162,17 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, push!(args, options[:δ]) # termination function - if invariant == nothing - termination = (k, set, t0) -> termination_N(N, k, set, t0) - else - termination = - (k, set, t0) -> termination_inv_N(N, invariant, k, set, t0) - end + termination = get_termination_function(N, invariant) push!(args, termination) + info("- Computing successors") + + # progress meter + progress_meter = (N != nothing) ? + Progress(N, 1, "Computing successors ") : + nothing + push!(args, progress_meter) + # choose algorithm backend algorithm = options[:algorithm] if algorithm == "explicit" @@ -180,14 +185,16 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}}, push!(args, res) # call the adequate function with the given arguments list - info("- Computing successors") @timing begin index, skip = available_algorithms[algorithm_backend]["func"](args...) - if index < N || skip - # shrink result array - info("terminated prematurely, only computed $index/$N steps") - deleteat!(res, (skip ? index : index + 1):N) + end + + # shrink result array + if skip || index < N + if N != nothing + info("termination after only $index of $N steps") end + deleteat!(res, (skip ? index : index + 1):length(res)) end # return the result @@ -210,10 +217,41 @@ function reach(system::IVP{<:AbstractContinuousSystem}, return reach(Δ, invariant, options) end -function termination_N(N, k, set, t0) +function get_termination_function(N::Nothing, invariant::Nothing) + termination_function = (k, set, t0) -> termination_unconditioned() + warn("no termination condition specified; the reachability analysis will " * + "not terminate") + return termination_function +end + +function get_termination_function(N::Int, invariant::Nothing) + return (k, set, t0) -> termination_N(N, k, t0) +end + +function get_termination_function(N::Nothing, invariant::LazySet) + return (k, set, t0) -> termination_inv(invariant, set, t0) +end + +function get_termination_function(N::Int, invariant::LazySet) + return (k, set, t0) -> termination_inv_N(N, invariant, k, set, t0) +end + +function termination_unconditioned() + return (false, false) +end + +function termination_N(N, k, t0) return (k >= N, false) end +function termination_inv(inv, set, t0) + if isdisjoint(set, inv) + return (true, true) + else + return (false, false) + end +end + function termination_inv_N(N, inv, k, set, t0) if k >= N return (true, false) @@ -257,3 +295,6 @@ end function has_constant_directions(block_options, i::Int) return true end + +# no-op progress meter for unbounded time +function update!(::Nothing, ::Int) end diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach_blocks.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach_blocks.jl index b8a5b176..c1d1c4da 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach_blocks.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach_blocks.jl @@ -38,7 +38,8 @@ given block indices. @inline row(ϕpowerk::SparseMatrixExp, bi::Int) = Matrix(get_row(ϕpowerk, bi)) @inline block(ϕpowerk_πbi::AbstractMatrix, bj::UnitRange{Int}) = ϕpowerk_πbi[:, bj] @inline block(ϕpowerk_πbi::AbstractMatrix, bj::Int) = ϕpowerk_πbi[:, [bj]] -@inline store!(res, k, X, t0, t1, N) = (res[k] = ReachSet(X, t0, t1)) +@inline store!(res, k, X, t0, t1, N::Int) = (res[k] = ReachSet(X, t0, t1)) +@inline store!(res, k, X, t0, t1, N::Nothing) = (push!(res, ReachSet(X, t0, t1))) # sparse function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, @@ -47,21 +48,23 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, overapproximate::Function, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, output_function::Union{Function, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, δ::NUM, termination::Function, + progress_meter::Union{Progress, Nothing}, res::Vector{<:ReachSet} )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) array = CartesianProductArray(Xhat0[blocks]) X_store = (output_function == nothing) ? array : box_approximation(output_function(array)) t0 = zero(δ) t1 = δ - store!(res, 1, X_store, t0, t1, NUM) + store!(res, 1, X_store, t0, t1, N) terminate, skip = termination(1, X_store, t0) if terminate return 1, skip @@ -81,9 +84,8 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, end k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] Xhatk_bi = ZeroSet(length(bi)) @@ -104,7 +106,7 @@ function reach_blocks!(ϕ::SparseMatrixCSC{NUM, Int}, box_approximation(output_function(array)) t0 = t1 t1 += δ - store!(res, k, X_store, t0, t1, NUM) + store!(res, k, X_store, t0, t1, N) terminate, skip = termination(k, X_store, t0) if terminate @@ -133,21 +135,23 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, overapproximate::Function, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, output_function::Union{Function, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, δ::NUM, termination::Function, + progress_meter::Union{Progress, Nothing}, res::Vector{<:ReachSet} )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) array = CartesianProductArray(Xhat0[blocks]) X_store = (output_function == nothing) ? array : box_approximation(output_function(array)) t0 = zero(δ) t1 = δ - store!(res, 1, X_store, t0, t1, NUM) + store!(res, 1, X_store, t0, t1, N) terminate, skip = termination(1, X_store, t0) if terminate return 1, skip @@ -170,9 +174,8 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, arr_length = (U == nothing) ? length(partition) : length(partition) + 1 arr = Vector{LazySet{NUM}}(undef, arr_length) k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] for (j, bj) in enumerate(partition) @@ -192,7 +195,7 @@ function reach_blocks!(ϕ::AbstractMatrix{NUM}, box_approximation(output_function(array)) t0 = t1 t1 += δ - store!(res, k, X_store, t0, t1, NUM) + store!(res, k, X_store, t0, t1, N) terminate, skip = termination(k, X_store, t0) if terminate @@ -224,21 +227,23 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, overapproximate::Function, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, output_function::Union{Function, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, δ::NUM, termination::Function, + progress_meter::Union{Progress, Nothing}, res::Vector{<:ReachSet} )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) array = CartesianProductArray(Xhat0[blocks]) X_store = (output_function == nothing) ? array : box_approximation(output_function(array)) t0 = zero(δ) t1 = δ - store!(res, 1, X_store, t0, t1, NUM) + store!(res, 1, X_store, t0, t1, N) terminate, skip = termination(1, X_store, t0) if terminate return 1, skip @@ -258,9 +263,8 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, end k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] ϕpowerk_πbi = row(ϕpowerk, bi) @@ -286,7 +290,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, box_approximation(output_function(array)) t0 = t1 t1 += δ - store!(res, k, X_store, t0, t1, NUM) + store!(res, k, X_store, t0, t1, N) terminate, skip = termination(k, X_store, t0) if terminate @@ -309,21 +313,23 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, overapproximate::Function, overapproximate_inputs::Function, n::Int, - N::Int, + N::Union{Int, Nothing}, output_function::Union{Function, Nothing}, blocks::AbstractVector{Int}, partition::AbstractVector{<:Union{AbstractVector{Int}, Int}}, δ::NUM, termination::Function, + progress_meter::Union{Progress, Nothing}, res::Vector{<:ReachSet} )::Tuple{Int, Bool} where {NUM} + update!(progress_meter, 1) array = CartesianProductArray(Xhat0[blocks]) X_store = (output_function == nothing) ? array : box_approximation(output_function(array)) t0 = zero(δ) t1 = δ - store!(res, 1, X_store, t0, t1, NUM) + store!(res, 1, X_store, t0, t1, N) terminate, skip = termination(1, X_store, t0) if terminate return 1, skip @@ -345,9 +351,8 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, arr_length = (U == nothing) ? length(partition) : length(partition) + 1 arr = Vector{LazySet{NUM}}(undef, arr_length) k = 2 - p = Progress(N, 1, "Computing successors ") @inbounds while true - update!(p, k) + update!(progress_meter, k) for i in 1:b bi = partition[blocks[i]] ϕpowerk_πbi = row(ϕpowerk, bi) @@ -371,7 +376,7 @@ function reach_blocks!(ϕ::SparseMatrixExp{NUM}, box_approximation(output_function(array)) t0 = t1 t1 += δ - store!(res, k, X_store, t0, t1, NUM) + store!(res, k, X_store, t0, t1, N) terminate, skip = termination(k, X_store, t0) if terminate diff --git a/src/ReachSets/ReachSets.jl b/src/ReachSets/ReachSets.jl index b5b5e92c..10fc8ca7 100644 --- a/src/ReachSets/ReachSets.jl +++ b/src/ReachSets/ReachSets.jl @@ -10,7 +10,7 @@ using LazySets, MathematicalSystems, HybridSystems, Expokit, Optim, # fix namespace conflicts with MathematicalSystems using LazySets: LinearMap -using Reachability: info +using Reachability: info, warn include("../compat.jl") From 216fb4c536a682100cd7402a504ffca3e3218983 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 16 Mar 2019 17:33:03 +0100 Subject: [PATCH 2/3] fix name bug --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index 7c77632c..c1bbaf2a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -58,7 +58,7 @@ solve(system::AbstractSystem, options::Pair{Symbol,<:Any}...) = function solve!(problem::InitialValueProblem, options_input::Options; - op::ContinuousPost=default_operator(system), + op::ContinuousPost=default_operator(problem), invariant::Union{LazySet, Nothing}=nothing )::AbstractSolution From 29ab5464296be3fe98c814d3f466a7d60f5e8aae Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 16 Mar 2019 17:33:40 +0100 Subject: [PATCH 3/3] unit tests for unbounded time --- test/Reachability/solve_continuous.jl | 32 +++++++++++++++++---------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/test/Reachability/solve_continuous.jl b/test/Reachability/solve_continuous.jl index 12c099b1..5fac3fcc 100644 --- a/test/Reachability/solve_continuous.jl +++ b/test/Reachability/solve_continuous.jl @@ -121,15 +121,15 @@ s = solve(system, Options(:T=>0.1), # Affine ODE with a nondeterministic input, x' = Ax + Bu # ======================================================= # linear ODE: x' = Ax + Bu, u ∈ U -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 - 0.38318 0.616014 0.518412 0.778765] +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 + 0.38318 0.616014 0.518412 0.778765] B = [0.866688 0.771231 - 0.065935 0.349839 - 0.109643 0.904222 - 0.292474 0.227857] + 0.065935 0.349839 + 0.109643 0.904222 + 0.292474 0.227857] # non-deterministic inputs U = Interval(0.99, 1.01) × Interval(0.99, 1.01) @@ -142,10 +142,18 @@ U1 = ConstantInput(U) U2 = U # use internal wrapping for inputs in [U1, U2] - # instantiate system - Δ = CLCCS(A, B, nothing, U) - 𝒮 = InitialValueProblem(Δ, X0) - # default options (computes all variables) - s = solve(𝒮, :T=>0.1) + s = solve(IVP(CLCCS(A, B, nothing, U), X0), :T=>0.1) end + +# ====================== +# Unbounded-time setting +# ====================== +X = LinearConstraint([1., 1., 0., 0.], 9.5) +property = SafeStatesProperty(LinearConstraint([1., 1., 0., 0.], 10.)) + +# default options (computes all variables) +s = Reachability.solve!(IVP(CLCCS(A, B, X, U), X0), Options(:T=>Inf); invariant=X) # temporary workaround for invariant (#507) +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)