Skip to content

Commit

Permalink
Merge pull request #143 from JuliaReach/mforets/algos_input
Browse files Browse the repository at this point in the history 
Update algorithms for reachability with inputs
  • Loading branch information
@mforets
mforets authored Apr 30, 2020
2 parents dbb9c4b + 04b7e7d commit 548d097
Show file tree
Hide file tree
Showing 13 changed files with 336 additions and 41 deletions.
2 changes: 1 addition & 1 deletion src/Algorithms/ASB07/post.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ function post(alg::ASB07, ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs...
U = inputset(ivp_discr)
@assert isa(U, LazySet) "expcted input of type `<:LazySet`, but got $(typeof(U))"
U = _convert_or_overapproximate(Zonotope, U)
reach_inhomog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, reduction_method)
reach_inhomog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, recursive, reduction_method)
end

return Flowpipe(F)
Expand Down
39 changes: 38 additions & 1 deletion src/Algorithms/ASB07/reach_homog.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
return F
end

# non-recursiv implementation; to get more accurate interval matrix powers Φ^k
# non-recursive implementation; to get more accurate interval matrix powers Φ^k
# we use the IntervalMatrices.IntervalMatrixPower interface
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
Ω0::Zonotope{N, VN, MN},
Expand Down Expand Up @@ -68,3 +68,40 @@ function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
end
return F
end

# case with an invariant
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
Ω0::Zonotope{N, VN, MN},
Φ::AbstractMatrix,
NSTEPS::Integer,
δ::N,
max_order::Integer,
X::LazySet,
recursive::Val{true},
reduction_method::AbstractReductionMethod) where {N, VN, MN}
# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# split the interval matrix into center and radius
Φc, Φs = _split(Φ)

k = 2
@inbounds while k <= NSTEPS
Zₖ₋₁ = set(F[k-1])
cₖ₋₁ = Zₖ₋₁.center
Gₖ₋₁ = Zₖ₋₁.generators

Zₖ = _overapproximate_interval_linear_map(Φc, Φs, cₖ₋₁, Gₖ₋₁)
Zₖʳ = _reduce_order(Zₖ, max_order, reduction_method)
_is_intersection_empty(X, Zₖʳ) && break

k += 1
Δt += δ
F[k] = ReachSet(Zₖʳ, Δt)
end
if k < NSTEPS + 1
resize!(F, k-1)
end
return F
end
85 changes: 85 additions & 0 deletions src/Algorithms/ASB07/reach_inhomog.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
# case with input and without invariant
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
Ω0::Zonotope{N, VN, MN},
Φ::AbstractMatrix,
NSTEPS::Integer,
δ::N,
max_order::Integer,
X::Universe,
U::Zonotope,
recursive::Val{true},
reduction_method::AbstractReductionMethod) where {N, VN, MN}
# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# input sequence
Wk₊ = copy(U)

# split the interval matrix into center and radius
Φc, Φs = _split(Φ)

k = 2
@inbounds while k <= NSTEPS
Zₖ₋₁ = set(F[k-1])
cₖ₋₁ = Zₖ₋₁.center
Gₖ₋₁ = Zₖ₋₁.generators

Zₖ = _overapproximate_interval_linear_map(Φc, Φs, cₖ₋₁, Gₖ₋₁)
Zₖ = _minkowski_sum(Wk₊, Zₖ)
Zₖʳ = _reduce_order(Zₖ, max_order, reduction_method)

k += 1
Δt += δ
F[k] = ReachSet(Zₖʳ, Δt)

Wk₊ = _overapproximate_interval_linear_map(Φc, Φs, Wk₊.center, Wk₊.generators)
Wk₊ = _reduce_order(Wk₊, max_order, reduction_method)
end
return F
end

# case with input and with invariant
function reach_homog_ASB07!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
Ω0::Zonotope{N, VN, MN},
Φ::AbstractMatrix,
NSTEPS::Integer,
δ::N,
max_order::Integer,
X::LazySet,
U::Zonotope,
recursive::Val{true},
reduction_method::AbstractReductionMethod) where {N, VN, MN}
# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# input sequence
Wk₊ = copy(U)

# split the interval matrix into center and radius
Φc, Φs = _split(Φ)

k = 2
@inbounds while k <= NSTEPS
Zₖ₋₁ = set(F[k-1])
cₖ₋₁ = Zₖ₋₁.center
Gₖ₋₁ = Zₖ₋₁.generators

Zₖ = _overapproximate_interval_linear_map(Φc, Φs, cₖ₋₁, Gₖ₋₁)
Zₖ = _minkowski_sum(Wk₊, Zₖ)
Zₖʳ = _reduce_order(Zₖ, max_order, reduction_method)
_is_intersection_empty(X, Zₖʳ) && break

k += 1
Δt += δ
F[k] = ReachSet(Zₖʳ, Δt)

Wk₊ = _overapproximate_interval_linear_map(Φc, Φs, Wk₊.center, Wk₊.generators)
Wk₊ = _reduce_order(Wk₊, max_order, reduction_method)
end
if k < NSTEPS + 1
resize!(F, k-1)
end
return F
end
39 changes: 23 additions & 16 deletions src/Algorithms/BOX/reach_homog.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,7 @@ function reach_homog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
δ::N,
X::Universe,
recursive::Val{false}) where {N, VNC, VNR}
#=
println("initial center: $(Ω0.center)")
println("initial radius: $(Ω0.radius)")
println("matrix: $Φ")
println("NSTEPS = $NSTEPS")
=#

# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)
Expand Down Expand Up @@ -95,22 +90,34 @@ function reach_homog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# preallocations
c = Vector{VNC}(undef, NSTEPS)
r = Vector{VNR}(undef, NSTEPS)

c[1] = Ω0.center
r[1] = Ω0.radius

# cache for powers of Φ
Φ_power_k = _get_cache(Φ)
Φ_power_k_cache = similar(Φ_power_k)

k = 2
@inbounds while k <= NSTEPS
Rₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
Δt += δ
F[k] = ReachSet(Rₖ, Δt)

mul!(Φ_power_k_cache, Φ_power_k, Φ)
copyto!(Φ_power_k, Φ_power_k_cache)
k += 1
# Hₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
c[k] = Φ_power_k * c[1]
r[k] = abs.(Φ_power_k) * r[1]
Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

_is_intersection_empty(X, Hₖ) && break
Δt += δ
F[k] = ReachSet(Hₖ, Δt)

mul!(Φ_power_k_cache, Φ_power_k, Φ)
copyto!(Φ_power_k, Φ_power_k_cache)
k += 1
end
if k < NSTEPS + 1
resize!(F, k-1)
end
return F
end

# homogeneous case with an invariant, recursive implementation
# not implemented
99 changes: 96 additions & 3 deletions src/Algorithms/BOX/reach_inhomog.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,102 @@
# ===================

# inhomogeneous case without invariant, non-recursive implementation
function reach_inhomog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
Ω0::Hyperrectangle{N, VNC, VNR},
Φ::AbstractMatrix,
NSTEPS::Integer,
δ::N,
X::Universe,
U::Hyperrectangle,
recursive::Val{false}) where {N, VNC, VNR}

# inhomogeneous case without invariant, recursive implementation
# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# inhomogeneous case with an invariant, non-recursive implmentation
# preallocations
c = Vector{VNC}(undef, NSTEPS)
r = Vector{VNR}(undef, NSTEPS)

# inhomogeneous case with an invariant, recursive implementation
c[1] = Ω0.center
r[1] = Ω0.radius

# input sequence
Wk₊ = copy(U)

# cache for powers of Φ
Φ_power_k = _get_cache(Φ)
Φ_power_k_cache = similar(Φ_power_k)

k = 2
@inbounds while k <= NSTEPS
# Hₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
c[k] = Φ_power_k * c[1] + Wk₊.center
r[k] = abs.(Φ_power_k) * r[1] + Wk₊.radius
Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

Δt += δ
F[k] = ReachSet(Hₖ, Δt)

Wk₊.center .+= Φ_power_k * U.center
Wk₊.radius .+= Φ_power_k * U.radius

mul!(Φ_power_k_cache, Φ_power_k, Φ)
copyto!(Φ_power_k, Φ_power_k_cache)
k += 1
end
return F
end

# homogeneous case with an invariant, non-recursive implementation
function reach_inhomog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
Ω0::Hyperrectangle{N, VNC, VNR},
Φ::AbstractMatrix,
NSTEPS::Integer,
δ::N,
X::LazySet,
U::Hyperrectangle,
recursive::Val{false}) where {N, VNC, VNR}

# initial reach set
Δt = zero(N) .. δ
@inbounds F[1] = ReachSet(Ω0, Δt)

# preallocations
c = Vector{VNC}(undef, NSTEPS)
r = Vector{VNR}(undef, NSTEPS)

c[1] = Ω0.center
r[1] = Ω0.radius

# input sequence
Wk₊ = copy(U)

# cache for powers of Φ
Φ_power_k = _get_cache(Φ)
Φ_power_k_cache = similar(Φ_power_k)

k = 2
@inbounds while k <= NSTEPS
# Hₖ = overapproximate(Φ_power_k * Ω0, Hyperrectangle)
c[k] = Φ_power_k * c[1] + Wk₊.center
r[k] = abs.(Φ_power_k) * r[1] + Wk₊.radius
Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

_is_intersection_empty(X, Hₖ) && break

Δt += δ
F[k] = ReachSet(Hₖ, Δt)

Wk₊.center .+= Φ_power_k * U.center
Wk₊.radius .+= Φ_power_k * U.radius

mul!(Φ_power_k_cache, Φ_power_k, Φ)
copyto!(Φ_power_k, Φ_power_k_cache)
k += 1
end
if k < NSTEPS + 1
resize!(F, k-1)
end
return F
end
2 changes: 1 addition & 1 deletion src/Algorithms/GLGM06/GLGM06.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ end

# convenience constructor using symbols
function GLGM06(; δ::N,
approx_model::AM=CorrectionHull(order=10, exp=:base),
approx_model::AM=Forward(sih=:concrete, exp=:base, setops=:lazy), #CorrectionHull(order=10, exp=:base),
max_order::Int=5,
static::Bool=false,
dim::Union{Int, Missing}=missing,
Expand Down
1 change: 1 addition & 0 deletions src/Algorithms/GLGM06/post.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ function post(alg::GLGM06, ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs..

reach_homog_GLGM06!(F, Ω0, Φ, NSTEPS, δ, max_order, X, preallocate)
else
# TODO: implement preallocate option for this scenario
U = inputset(ivp_discr)
@assert isa(U, LazySet) "expcted input of type `<:LazySet`, but got $(typeof(U))"
U = _convert_or_overapproximate(Zonotope, U)
Expand Down
Loading

0 comments on commit 548d097

Please sign in to comment.