Merge pull request #340 from JuliaReach/mforets/update_basics

Add exponential LGG09 w/ time-varying input sets

docs/src/man/basics.md

@@ -50,14 +50,58 @@ discrete modes that are compatible with the dynamics.
 
 ## Set representations
 
-### Convex sets
+### Representing sets with support functions
+
+Support functions are one of the central tools in set-based reachability. The *support function* of a closed and bounded convex set $\mathcal{X} \subseteq \mathbb{R}^n$ attributes to a direction vector $d \in \mathbb{R}^n$ the real number
+
+```math
+  ρ(d, \mathcal{X}) = \max \{ d^T x : x \in \mathcal{X} \},
+```
+where $d^T$ denotes the transpose of the (column) vector $d$.
+
+TODO: example with triangle
+A convex set can be represe
+
+The
 
 ### Polyhedra
 
 ### Zonotopes
 
+Zonotopes are a sub-class of polytopes defined as the image of a unit cube under
+
+
+### Hausdorff distance
+
+
+
 ### LazySets
 
+
+### Taylor models
+
+## Taylor models
+
+
+```@example
+
+using ReachabilityAnalysis, Plots
+
+f(x) = -6x^3 + (13/3)x^2 + (31/3)x
+dom = -3.5 .. 3.5
+
+plot(f, -3.5, 3.5, lab="f", xlab="x")
+
+x = Taylor1(5)
+set_taylor1_varname("x")
+f(x)
+
+rem = 0 .. 0
+x0 = 0.0
+dom = -3.5 .. 3.5
+tm = TaylorModel1(f(x), rem, x0, dom)
+```
+
 ## Reach-sets
 
 We consider as a running example in this section the simple harmonic oscillator,

src/Algorithms/LGG09/reach_homog.jl

@@ -115,70 +115,52 @@ end
 # Functionality that requires ExponentialUtilities.jl
 # ------------------------------------------------------------
 
-function load_exponential_utilities_LGG09()
+function load_krylov_LGG09_homog()
 return quote
 
-# recursive version, default expv
-function reach_homog_dir_LGG09_expv!(out, Ω₀, Aᵀ, ℓ, NSTEPS, recursive::Val{:true})
-    rᵢ = copy(ℓ)
-    rᵢ₊₁ = similar(rᵢ)
-
-    @inbounds for i in 1:NSTEPS
-        out[i] = ρ(rᵢ, Ω₀)
-
-        # update cache for the next iteration
-        rᵢ₊₁ = expv(1.0, Aᵀ, rᵢ) # computes exp(Aᵀ * 1.0) * rᵢ
-        copy!(rᵢ, rᵢ₊₁)
-    end
-    return out
-end
-
-# ρ(ℓ, Ω₀), ρ(exp(Aᵀ) * ℓ, Ω₀), ρ(exp(2Aᵀ) * ℓ, Ω₀)
-function reach_homog_dir_LGG09_expv!(out, Ω₀, Aᵀ, ℓ, NSTEPS, recursive::Val{:false})
-    rᵢ = deepcopy(ℓ) # if ℓ is a sev => this is a sev
-
-    @inbounds for i in 1:NSTEPS
-        out[i] = ρ(rᵢ, Ω₀)
-
-        # update cache for the next iteration
-        rᵢ = expv(i*1.0, Aᵀ, ℓ)
-    end
-    return out
-end
-
-# non-recursive version using precomputed Krylov subspace
-function reach_homog_dir_LGG09_expv_pk!(out, Ω₀, Aᵀ, ℓ, NSTEPS, recursive::Val{:false})
-    rᵢ = deepcopy(ℓ) # if ℓ is a sev => this is a sev
-    Ks = arnoldi(Aᵀ, ℓ, tol=1e-18)
-
-    @inbounds for i in 1:NSTEPS
-        out[i] = ρ(rᵢ, Ω₀)
-
-        # update cache for the next iteration
-        expv!(rᵢ, i*1.0, Ks)
-    end
-    return out
-end
-
-# this function computes the sequence
-# ``ρ(ℓ, Ω₀)``, ``ρ(exp(Aᵀ) * ℓ, Ω₀)``, ``ρ(exp(2Aᵀ) * ℓ, Ω₀)`` until ``ρ(exp(NSTEPS * Aᵀ) * ℓ, Ω₀)``
-function reach_homog_dir_LGG09_expv_pk2!(out, Ω₀, Aᵀ, ℓ, NSTEPS, recursive::Val{:false};
-                                        hermitian=false, m=min(30, size(Aᵀ, 1)), tol=1e-7)
-
-    TA, Tb = eltype(Aᵀ), eltype(ℓ)
+# Compute the sequence:
+#
+# ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀), ρ(ℓ, Φ^2 Ω₀), ρ(ℓ, Φ^3 Ω₀), ...
+#
+# Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
+# the direction ℓ.
+#
+# Method (see [1]):
+#
+# out[1] <- ρ(ℓ, Ω₀)
+#
+# out[2] <- ρ(ℓ, Φ Ω₀) = ρ(Φᵀ ℓ, Ω₀)
+#
+# out[3] <- ρ(ℓ, Φ^2 Ω₀) = ρ((Φᵀ)^2 ℓ, Ω₀)
+#
+# out[4] <- ρ(ℓ, Φ^3 Ω₀) = ρ((Φᵀ)^3 ℓ, Ω₀)
+#
+# and so on.
+#
+# [1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
+#     High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
+#     Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
+#     Proceedings of the 21st International Conference on Hybrid Systems: Computation
+#     and Control: 41–50.
+#
+function reach_homog_krylov_LGG09!(out, Ω₀::LazySet, Aᵀδ::AbstractMatrix,
+                                   ℓ::AbstractVector, NSTEPS;
+                                   hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)
+
+    # initialization of the krylov subspace
+    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
     T = promote_type(TA, Tb)
     Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
-    arnoldi!(Ks, Aᵀ, ℓ; m=m, ishermitian=hermitian, tol=tol)
+    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)
 
+    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
     rᵢ = deepcopy(ℓ)
 
     @inbounds for i in 1:NSTEPS
         out[i] = ρ(rᵢ, Ω₀)
-
-        # update cache for the next iteration
         expv!(rᵢ, i*1.0, Ks)
     end
     return out
 end
 
-end end  # quote / load_exponential_utilities_LGG09()
+end end  # quote / load_krylov_LGG09_homog()

src/Algorithms/LGG09/reach_inhomog.jl

@@ -276,3 +276,113 @@ function reach_inhomog_LGG09!(F::Vector{RT},
     return ρℓ
 end
 =#
+
+# ------------------------------------------------------------
+# Functionality that requires ExponentialUtilities.jl
+# ------------------------------------------------------------
+
+function load_krylov_LGG09_inhomog()
+return quote
+
+# Compute the sequence with constant input sets:
+#
+# ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀ ⊕ V), ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V ⊕ V), ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V ⊕ Φ V ⊕ V), ...
+#
+# Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
+# the direction ℓ.
+#
+# Method (see [1]):
+#
+# out[1] <- ρ(ℓ, Ω₀)
+#
+# out[2] <- ρ(ℓ, Φ Ω₀ ⊕ V) = ρ(ℓ, Φ Ω0) + ρ(ℓ, V) = ρ(Φᵀ ℓ, Ω₀) + ρ(ℓ, V)
+#
+# out[3] <- ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V ⊕ V) = ρ((Φᵀ)^2 ℓ, Ω₀) + ρ(Φᵀ ℓ, V) + ρ(ℓ, V)
+#
+# out[4] <- ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V ⊕ Φ V ⊕ V) = ρ((Φᵀ)^3 ℓ, Ω₀) + ρ((Φᵀ)^2 ℓ, V) + ρ(Φᵀ ℓ, V) + ρ(ℓ, V)
+#
+# and so on.
+#
+# [1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
+#     High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
+#     Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
+#     Proceedings of the 21st International Conference on Hybrid Systems: Computation
+#     and Control: 41–50.
+#
+function reach_inhomog_krylov_LGG09!(out, Ω₀::LazySet, V::LazySet, Aᵀδ::AbstractMatrix,
+                                     ℓ::AbstractVector, NSTEPS;
+                                     hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)
+
+    # initialization of the krylov subspace
+    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
+    T = promote_type(TA, Tb)
+    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
+    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)
+
+    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
+    rᵢ = deepcopy(ℓ)
+    out[1] = ρ(ℓ, Ω₀)
+
+    # accumulated support vector sum due to the inputs
+    s = zero(T)
+
+    @inbounds for i in 1:NSTEPS-1
+        s += ρ(rᵢ, V)
+        expv!(rᵢ, i*1.0, Ks)
+        out[i+1] = ρ(rᵢ, Ω₀) + s
+    end
+    return out
+end
+
+# Compute the sequence with time-varying input sets:
+#
+# ρ(ℓ, Ω₀), ρ(ℓ, Φ Ω₀ ⊕ V₀), ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V₀ ⊕ V₁), ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V₀ ⊕ Φ V₁ ⊕ V₂), ...
+#
+# Using Krylov subspace approximations to compute the action of Φ := exp(Aδ) over
+# the direction ℓ.
+#
+# Method (see [1]):
+#
+# out[1] <- ρ(ℓ, Ω₀)
+#
+# out[2] <- ρ(ℓ, Φ Ω₀ ⊕ V₀) = ρ(ℓ, Φ Ω0) + ρ(ℓ, V₀) = ρ(Φᵀ ℓ, Ω₀) + ρ(ℓ, V₀)
+#
+# out[3] <- ρ(ℓ, Φ^2 Ω₀ ⊕ Φ V₀ ⊕ V₁) = ρ((Φᵀ)^2 ℓ, Ω₀) + ρ(Φᵀ ℓ, V₀) + ρ(ℓ, V₁)
+#
+# out[4] <- ρ(ℓ, Φ^3 Ω₀ ⊕ Φ^2 V₀ ⊕ Φ V₁ ⊕ V₂) = ρ((Φᵀ)^3 ℓ, Ω₀) + ρ((Φᵀ)^2 ℓ, V₀) + ρ(Φᵀ ℓ, V₁) + ρ(ℓ, V₂)
+#
+# and so on.
+#
+# [1] Reach Set Approximation through Decomposition with Low-dimensional Sets and
+#     High-dimensional Matrices. Sergiy Bogomolov, Marcelo Forets, Goran Frehse,
+#     Frédéric Viry, Andreas Podelski and Christian Schilling (2018) HSCC'18
+#     Proceedings of the 21st International Conference on Hybrid Systems: Computation
+#     and Control: 41–50.
+#
+function reach_inhomog_krylov_LGG09!(out, Ω₀::LazySet, V::Vector{<:LazySet}, Aᵀδ::AbstractMatrix,
+                                     ℓ::AbstractVector, NSTEPS;
+                                     hermitian=false, m=min(30, size(Aᵀδ, 1)), tol=1e-7)
+
+    # initialization of the krylov subspace
+    TA, Tb = eltype(Aᵀδ), eltype(ℓ)
+    T = promote_type(TA, Tb)
+    Ks = KrylovSubspace{T, real(T)}(length(ℓ), m)
+    arnoldi!(Ks, Aᵀδ, ℓ; m=m, ishermitian=hermitian, tol=tol)
+
+    # rᵢ stores is the cache for each vector: (Φᵀ)^i ℓ
+    r = Vector{Vector{T}}(undef, NSTEPS)
+    r[1] = deepcopy(ℓ)
+    out[1] = ρ(ℓ, Ω₀)
+
+    @inbounds for i in 1:NSTEPS-1
+        s = zero(T)
+        for j in 1:i
+            s += ρ(r[i-j+1], V[j])
+        end
+        expv!(r[i+1], i*1.0, Ks)
+        out[i+1] = ρ(r[i+1], Ω₀) + s
+    end
+    return out
+end
+
+end end  # quote / load_krylov_LGG09_inhomog()

src/Flowpipes/flowpipes.jl

@@ -341,6 +341,7 @@ function project(fp::Flowpipe, vars::NTuple{D, T}) where {D, T<:Integer}
     end
 end
 
+project(fp::Flowpipe, vars::Int) = project(fp, (vars,))
 project(fp::Flowpipe, vars::AbstractVector{<:Int}) = project(fp, Tuple(vars))
 project(fp::Flowpipe; vars) = project(fp, Tuple(vars))
 project(fp::Flowpipe, i::Int, vars) = project(fp[i], vars)

src/Flowpipes/reachsets.jl

@@ -513,6 +513,11 @@ end
     vars_idx = indexin(variables, all_variables) |> Vector{Int}
 end
 
+# concrete projection onto a single variable
+function project(R::AbstractLazyReachSet, variable::Int; check_vars::Bool=true)
+    return project(R, (variable,), check_vars=check_vars)
+end
+
 # lazy projection of a reach-set
 function Projection(R::AbstractLazyReachSet, variables::NTuple{D, M},
                     check_vars::Bool=true) where {D, M<:Integer}

src/Initialization/init_ExponentialUtilities.jl

@@ -2,4 +2,5 @@ eval(quote
     using .ExponentialUtilities: expv, expv!, arnoldi, arnoldi!, KrylovSubspace
 end)
 
-eval(load_exponential_utilities_LGG09())
+eval(load_krylov_LGG09_homog())
+eval(load_krylov_LGG09_inhomog())