#507 - Remove invariant argument from post

src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl

@@ -251,17 +251,10 @@ function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂_input::Options)
     # convert matrix
     system = matrix_conversion(𝑆, 𝑂)
 
-    # apply transformation to the invariant
-    if hasmethod(stateset, Tuple{typeof(𝑆.s)})
-        invariant = 𝑆.s.X
-    else
-        invariant = Universe(𝑂_input[:n])
-    end
-
     if 𝑂[:mode] == "reach"
         info("Reachable States Computation...")
         @timing begin
-            Rsets = reach(𝑆, invariant, 𝑂)
+            Rsets = reach(𝑆, 𝑂)
             info("- Total")
         end
 

src/ReachSets/ContinuousPost/BFFPSV18/reach.jl

@@ -7,14 +7,13 @@ import LazySets.Approximations: overapproximate
 import ProgressMeter: update!
 
 """
-    reach(S, invariant, options)
+    reach(problem, options)
 
 Interface to reachability algorithms for an LTI system.
 
 ### Input
 
-- `S`         -- LTI system, discrete or continuous
-- `invariant` -- invariant
+- `problem`   -- initial value problem
 - `options`   -- additional options
 
 ### Output
@@ -28,27 +27,26 @@ A dictionary with available algorithms is available via
 
 The numeric type of the system's coefficients and the set of initial states
 is inferred from the first parameter of the system (resp. lazy set), ie.
-`NUM = first(typeof(S.s).parameters)`.
+`NUM = first(typeof(problem.s).parameters)`.
 """
-function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}},
-                        IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
-               invariant::Union{LazySet, Nothing},
+function reach(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
+                              IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
               )::Vector{<:ReachSet} where {NUM <: Real}
 
     # list containing the arguments passed to any reachability function
     args = []
 
     # coefficients matrix
-    A = S.s.A
+    A = problem.s.A
     push!(args, A)
 
     # determine analysis mode (sparse/dense) for lazy_expm mode
     if A isa SparseMatrixExp
         push!(args, Val(options[:assume_sparse]))
     end
 
-    n = statedim(S)
+    n = statedim(problem)
     blocks = options[:blocks]
     partition = convert_partition(options[:partition])
     T = options[:T]
@@ -58,11 +56,11 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}},
     if length(partition) == 1 && length(partition[1]) == n &&
             options[:block_options_init] == LinearMap
         info("- Skipping decomposition of X0")
-        Xhat0 = LazySet{NUM}[S.x0]
+        Xhat0 = LazySet{NUM}[problem.x0]
     else
         info("- Decomposing X0")
         @timing begin
-            Xhat0 = array(decompose(S.x0, options[:partition],
+            Xhat0 = array(decompose(problem.x0, options[:partition],
                                     options[:block_options_init]))
         end
     end
@@ -90,8 +88,8 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}},
     push!(args, Xhat0)
 
     # inputs
-    if !options[:assume_homogeneous] && inputdim(S) > 0
-        U = inputset(S)
+    if !options[:assume_homogeneous] && inputdim(problem) > 0
+        U = inputset(problem)
     else
         U = nothing
     end
@@ -162,6 +160,7 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}},
     push!(args, options[:δ])
 
     # termination function
+    invariant = stateset(problem.s)
     termination = get_termination_function(N, invariant)
     push!(args, termination)
 
@@ -201,30 +200,28 @@ function reach(S::Union{IVP{<:LDS{NUM}, <:LazySet{NUM}},
     return res
 end
 
-function reach(system::IVP{<:AbstractContinuousSystem},
-               invariant::Union{LazySet, Nothing},
+function reach(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}},
+                              IVP{<:CLCCS{NUM}, <:LazySet{NUM}}},
                options::TwoLayerOptions
-              )::Vector{<:ReachSet}
-    # ===================
-    # Time discretization
-    # ===================
+              )::Vector{<:ReachSet} where {NUM <: Real}
+
     info("Time discretization...")
-    Δ = @timing discretize(system, options[:δ], algorithm=options[:discretization],
-                                                exp_method=options[:exp_method],
-                                                sih_method=options[:sih_method])
+    Δ = @timing discretize(problem, options[:δ], algorithm=options[:discretization],
+                                                 exp_method=options[:exp_method],
+                                                 sih_method=options[:sih_method])
 
     Δ = matrix_conversion(Δ, options)
-    return reach(Δ, invariant, options)
+    return reach(Δ, options)
 end
 
-function get_termination_function(N::Nothing, invariant::Nothing)
+function get_termination_function(N::Nothing, invariant::Universe)
     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)
+function get_termination_function(N::Int, invariant::Universe)
     return (k, set, t0) -> termination_N(N, k, t0)
 end
 

src/ReachSets/discretize.jl

@@ -1,8 +1,3 @@
-const LDS = LinearDiscreteSystem
-const CLCDS = ConstrainedLinearControlDiscreteSystem
-
-@inline Id(n) = Matrix(1.0I, n, n)
-
 """
     discretize(𝑆, δ; [algorithm], [exp_method], [sih_method], [set_operations])
 
@@ -416,7 +411,7 @@ function discretize_firstorder(𝑆::InitialValueProblem,
         α = κ * RX0
         □α = Ballp(p, zeros(n), α)
         Ω0 = ConvexHull(X0, ϕ * X0 ⊕ □α)
-        return IVP(LDS(ϕ), Ω0)
+        return IVP(CLDS(ϕ, stateset(𝑆.s)), Ω0)
 
     elseif isaffine(𝑆)
         Uset = inputset(𝑆)
@@ -435,7 +430,7 @@ function discretize_firstorder(𝑆::InitialValueProblem,
             # transformation of the inputs
             □β = Ballp(p, zeros(n), β)
             Ud = ConstantInput(δ*U ⊕ □β)
-            return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0)
+            return IVP(CLCDS(ϕ, I(n), stateset(𝑆.s), Ud), Ω0)
 
         elseif Uset isa VaryingInput
             Ud = Vector{LazySet}(undef, length(Uset))
@@ -457,7 +452,7 @@ function discretize_firstorder(𝑆::InitialValueProblem,
                 Ud[i] = δ*Ui ⊕ □β
             end
             Ud = VaryingInput(Ud)
-            return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0)
+            return IVP(CLCDS(ϕ, I(n), stateset(𝑆.s), Ud), Ω0)
         else
             throw(ArgumentError("input of type $(typeof(U)) is not allowed"))
         end
@@ -525,6 +520,7 @@ function  discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSy
 
     # unwrap coefficient matrix and initial states
     A, X0 = 𝑆.s.A, 𝑆.x0
+    n = size(A, 1)
 
     # compute matrix ϕ = exp(Aδ)
     ϕ = exp_Aδ(A, δ, exp_method=exp_method)
@@ -534,15 +530,15 @@ function  discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSy
 
     # early return for homogeneous systems
     if inputdim(𝑆) == 0
-        return IVP(LDS(ϕ), Ω0)
+        return IVP(CLDS(ϕ, stateset(𝑆.s)), Ω0)
     end
 
     # compute matrix to transform the inputs
     Phi1Adelta = Φ₁(A, δ, exp_method=exp_method)
     U = inputset(𝑆)
     Ud = map(ui -> Phi1Adelta * ui, U)
 
-    return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0)
+    return IVP(CLCDS(ϕ, I(n), stateset(𝑆.s), Ud), Ω0)
 end
 
 """
@@ -646,7 +642,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous
     # early return for homogeneous systems
     if inputdim(𝑆) == 0
         Ω0 = _discretize_interpolation_homog(X0, ϕ, Einit, Val(set_operations))
-        return IVP(LDS(ϕ), Ω0)
+        return IVP(CLDS(ϕ, stateset(𝑆.s)), Ω0)
     end
 
     U = inputset(𝑆)
@@ -655,7 +651,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous
     Eψ0 = sih(Phi2Aabs * sih(A * U0))
     Ω0, Ud = _discretize_interpolation_inhomog(δ, U0, U, X0, ϕ, Einit, Eψ0, A, sih, Phi2Aabs, Val(set_operations))
 
-    return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0)
+    return IVP(CLCDS(ϕ, I(size(A, 1)), stateset(𝑆.s), Ud), Ω0)
 end
 
 # version using lazy sets and operations

src/Reachability.jl

@@ -11,12 +11,6 @@ using Reexport, RecipesBase, Memento, MathematicalSystems, HybridSystems,
 import LazySets.use_precise_ρ
 import LazySets.LinearMap
 
- # common aliases for MathematicalSystems types
- const CLCCS = ConstrainedLinearControlContinuousSystem
- const CLCDS = ConstrainedLinearControlDiscreteSystem
- const LCS = LinearContinuousSystem
- const LDS = LinearDiscreteSystem
-
 include("logging.jl")
 include("Options/dictionary.jl")
 include("Options/validation.jl")

src/Utils/normalization.jl

@@ -1,12 +1,3 @@
-# linear ivp in canonical form x' = Ax
-const LCF = LCS{N, AN} where {N, AN<:AbstractMatrix{N}}
-
-# affine ivp with constraints in canonical form x' = Ax + u, u ∈ U, x ∈ X
-const ACF = CLCCS{N, AN, IdentityMultiple{N}, XT, UT} where {N, AN<:AbstractMatrix{N}, XT<:LazySet{N}, UT<:AbstractInput}
-
-# type union of canonical forms
-const CF = Union{<:LCF, <:ACF}
-
 # accepted types of non-deterministic inputs (non-canonical form)
 const UNCF = Union{<:LazySet{N}, Vector{<:LazySet{N}}, <:AbstractInput} where {N}
 
@@ -30,21 +21,22 @@ system otherwise.
 ### Notes
 
 The normalization procedure consists of transforming a given system type into a
-"canonical" format that is used internally. The type union `CF` defines the
-systems which are considered canonical, i.e. which do not require normalization.
-More details are given below.
+"canonical" format that is used internally. More details are given below.
 
 ### Algorithm
 
 The implementation of `normalize` exploits `MathematicalSystems`'s' types, which
 carry information about the problem as a type parameter.
 
-Homogeneous ODEs of the form ``x' = Ax`` are canonical if the associated
-problem is a `LinearContinuousSystem` and `A` is a matrix.
+Homogeneous ODEs of the form ``x' = Ax, x ∈ \\mathcal{X}`` are canonical if the
+associated problem is a `ConstrainedLinearContinuousSystem` and `A` is a matrix.
+This type does not handle non-deterministic inputs.
+
 Note that a `LinearContinuousSystem` does not consider constraints on the
 state-space (such as an invariant); to specify state constraints, use a
-`ConstrainedLinearControlContinuousSystem`. Moreover, this type does not handle
-non-deterministic inputs.
+`ConstrainedLinearContinuousSystem`. If the passed system is a `LinearContinuousSystem`
+(i.e. no constraints) then the normalization fixes a universal set (`Universe`)
+as the constraint set.
 
 The generalization to canonical systems with constraints and possibly time-varying
 non-deterministic inputs is considered next. These systems are of the form
@@ -69,31 +61,72 @@ function normalize(system::AbstractSystem)
     throw(ArgumentError("the system type $(typeof(system)) is currently not supported"))
 end
 
-# systems already in canonical form, i.e. which don't need normalization
-normalize(system::CF) = system
+# x' = Ax, in the continuous case
+# x+ = Ax, in the discrete case
+function normalize(system::LCS{N, AN}) where {N, AN<:AbstractMatrix{N}}
+    n = statedim(system)
+    X = Universe(n)
+    return CLCS(system.A, X)
+end
+
+function normalize(system::LDS{N, AN}) where {N, AN<:AbstractMatrix{N}}
+    n = statedim(system)
+    X = Universe(n)
+    return CLDS(system.A, X)
+end
+
+# x' = Ax, x ∈ X in the continuous case
+# x+ = Ax, x ∈ X in the discrete case
+for CL_S in (:CLCS, :CLDS)
+    @eval begin
+        function normalize(system::$CL_S{N, AN, ST}) where {N, AN<:AbstractMatrix{N}, ST<:LazySet{N}}
+            n = statedim(system)
+            X = _wrap_invariant(stateset(system), n)
+            return $CL_S(system.A, X)
+        end
+    end
+end
+
+# x' = Ax + u, x ∈ X, u ∈ U in the continuous case
+# x+ = Ax + u, x ∈ X, u ∈ U in the discrete case
+for CLC_X in (:CLCCS, :CLCDS)
+    @eval begin
+        function normalize(system::$CLC_X{N, IdentityMultiple{N}, ST}) where {N, ST<:LazySet{N}}
+            λ = system.B.M.λ
+            if λ == oneunit(λ)
+                return system
+            else
+                n = statedim(system)
+                X = _wrap_invariant(stateset(system), n)
+                U = _wrap_inputs(system.U, Diagonal(fill(λ, n)))
+                return $CLC_X(system.A, I(n, N), X, U)
+            end
+        end
+    end
+end
 
 # x' = Ax + Bu, x ∈ X, u ∈ U in the continuous case
 # x+ = Ax + Bu, x ∈ X, u ∈ U in the discrete case
-for CLCS in (:CLCCS, :CLCDS)
+for CLC_S in (:CLCCS, :CLCDS)
     @eval begin
-        function normalize(system::$CLCS{N, AN, BN, XT, UT}) where {N, AN<:AbstractMatrix{N}, BN<:AbstractMatrix{N}, XT<:XNCF, UT<:UNCF}
+        function normalize(system::$CLC_S{N, AN, BN, XT, UT}) where {N, AN<:AbstractMatrix{N}, BN<:AbstractMatrix{N}, XT<:XNCF, UT<:UNCF}
             n = statedim(system)
-            X = _wrap_invariant(system.X, n)
+            X = _wrap_invariant(stateset(system), n)
             U = _wrap_inputs(system.U, system.B)
-            $CLCS(system.A, I(n, N), X, U)
+            $CLC_S(system.A, I(n, N), X, U)
         end
     end
 end
 
 # x' = Ax + Bu + c, x ∈ X, u ∈ U in the continuous case
 # x+ = Ax + Bu + c, x ∈ X, u ∈ U in the discrete case
-for CACS in (:CACCS, :CACDS)
+for CAC_S in (:CACCS, :CACDS)
     @eval begin
-        function normalize(system::$CACS{N, AN, BN, VN, XT, UT}) where {N, AN<:AbstractMatrix{N}, BN<:AbstractMatrix{N}, VN<:AbstractVector{N}, XT<:XNCF, UT<:UNCF}
+        function normalize(system::$CAC_S{N, AN, BN, VN, XT, UT}) where {N, AN<:AbstractMatrix{N}, BN<:AbstractMatrix{N}, VN<:AbstractVector{N}, XT<:XNCF, UT<:UNCF}
             n = statedim(system)
             X = _wrap_invariant(system.X, n)
             U = _wrap_inputs(system.U, system.B, system.c)
-            $CACS(system.A, I(n, N), X, U)
+            $CAC_S(system.A, I(n, N), X, U)
         end
     end
 end