#507 - Remove invariant argument from post

src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl

@@ -235,23 +235,29 @@ function init!(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)
 end
 
 """
-    post(𝒫::BFFPSV18, 𝑆::AbstractSystem, invariant, 𝑂::Options)
+    post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)
 
 Calculate the reachable states of the given initial value problem using `BFFPSV18`.
 
 ### Input
 
 - `𝒫` -- post operator of type `BFFPSV18`
 - `𝑆` -- sytem, initial value problem for a continuous ODE
-- `invariant` -- constraint invariant on the mode
 - `𝑂` -- algorithm-specific options
 """
-function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, invariant, 𝑂_input::Options)
+function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂_input::Options)
     𝑂 = TwoLayerOptions(merge(𝑂_input, 𝒫.options.specified), 𝒫.options.defaults)
 
     # 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

src/Reachability.jl

@@ -18,21 +18,18 @@ import LazySets.LinearMap
  const LDS = LinearDiscreteSystem
 
 include("logging.jl")
-include("Utils/Utils.jl")
 include("Options/dictionary.jl")
 include("Options/validation.jl")
 include("Options/default_options.jl")
+include("Utils/Utils.jl")
 include("Properties/Properties.jl")
 include("ReachSets/ReachSets.jl")
-include("Transformations/Transformations.jl")
 
 @reexport using Reachability.Utils,
                 Reachability.ReachSets,
-                Reachability.Properties,
-                Reachability.Transformations
+                Reachability.Properties
 
-export project,
-       Transformations
+export project
 
 include("solve.jl")
 include("plot_recipes.jl")

src/Utils/Utils.jl

@@ -8,7 +8,7 @@ using LazySets, MathematicalSystems, HybridSystems
 
 include("../compat.jl")
 
-import Reachability.@timing
+using Reachability: Options
 
 # Visualization
 export print_sparsity,
@@ -41,6 +41,9 @@ include("systems.jl")
 # normalization
 include("normalization.jl")
 
+# coordinate transformations
+include("transformations.jl")
+
 # abstract solution type
 include("AbstractSolution.jl")
 

src/Utils/transformations.jl

@@ -0,0 +1,84 @@
+export transform
+
+"""
+    transform(problem::InitialValueProblem, options::Options)
+
+Interface function that calls the respective transformation function.
+
+### Input
+
+- `problem` -- discrete or continuous initial-value problem
+- `option`  -- problem options
+
+### Output
+
+A tuple containing the transformed problem and the transformed options.
+
+### Notes
+
+The functions that are called in the background should return a the transformed
+system components `A`, `X0`, and `U`, and also an inverse transformation matrix `M`.
+If the system has an invariant, it is transformed as well.
+"""
+function transform(problem::InitialValueProblem, options::Options)
+    method = options[:coordinate_transformation]
+    if method == ""
+        nothing # no-op
+    elseif method == "schur"
+        options[:transformation_matrix] = Tm
+        problem = schur_transform(problem)
+    else
+        error("the transformation method $method is undefined")
+    end
+    return (problem, options)
+end
+
+"""
+    schur_transform(problem::InitialValueProblem)
+
+Applies a Schur transformation to a discrete or continuous initial-value problem.
+
+### Input
+
+- `S` -- discrete or continuous initial-value problem
+
+### Output
+
+Transformed problem.
+
+### Algorithm
+
+We use Julia's default `schurfact` function to compute a
+[Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)
+of the coefficients matrix ``A``.
+"""
+function schur_transform(problem::InitialValueProblem{PT, ST}
+                         ) where {PT <: Union{ConstrainedLinearControlDiscreteSystem, ConstrainedLinearControlContinuousSystem}, ST<:LazySet}
+
+    n = size(problem.s.A, 1)
+
+    # full (dense) matrix is required by schur
+    A_new, T_new = schur(Matrix(problem.s.A))
+
+    # recall that for Schur matrices, inv(T) == T'
+    Z_inverse = copy(transpose(T_new))
+
+    # apply transformation to the initial states
+    X0_new = Z_inverse * problem.x0
+
+    # apply transformation to the inputs
+    U_new = Z_inverse * problem.s.U
+
+    # obtain the transformation matrix for reverting the transformation again
+    T_inverse = F[:vectors]
+
+    # apply transformation to the invariant
+    if hasmethod(stateset, Tuple{typeof(problem.s)})
+        invariant_new = T_inverse * problem.s.X
+    else
+        invariant_new = Universe(n)
+    end
+
+    system_new = (A_new, I(n), invariant_new, U_new)
+    return InitialValueProblem(PT(system_new), X0_new)
+end

src/solve.jl

@@ -57,25 +57,21 @@ solve(system::AbstractSystem, options::Pair{Symbol,<:Any}...) =
     solve(system, Options(Dict{Symbol,Any}(options)))
 
 function solve!(problem::InitialValueProblem,
-                options_input::Options;
-                op::ContinuousPost=default_operator(problem),
-                invariant::Union{LazySet, Nothing}=nothing
+                options::Options;
+                op::ContinuousPost=default_operator(problem)
                )::AbstractSolution
 
-    system = normalize(problem.s)
-    problem = IVP(system, problem.x0)
-    options = init!(op, problem, options_input)
-
-    # coordinate transformation
-    if options[:coordinate_transformation] != ""
-        info("Transformation...")
-        (problem, transformation_matrix) =
-            @timing transform(problem, options[:coordinate_transformation])
-        options[:transformation_matrix] = transformation_matrix
-        invariant = options[:coordinate_transformation] * invariant
-    end
+    # normalize system to a canonical form if needed
+    problem = IVP(normalize(problem.s), problem.x0)
+
+    # initialize the algorithm-specific options
+    options = init!(op, problem, options)
+
+    # compute a coordinate transformation if needed
+    problem, options = transform(problem, options)
 
-    post(op, problem, invariant, options)
+    # run the continuous post-operator
+    post(op, problem, options)
 end
 
 """
@@ -167,8 +163,7 @@ function solve!(system::InitialValueProblem{<:HybridSystem,
         end
         reach_tube = solve!(IVP(loc, X0.X),
                             options_copy,
-                            op=opC,
-                            invariant=source_invariant)
+                            op=opC)
         inout_map = reach_tube.options[:inout_map]  # TODO temporary hack
         # get the property for the current location
         property_loc = property isa Dict ?

test/Reachability/solve_continuous.jl

@@ -153,7 +153,7 @@ 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 = Reachability.solve!(IVP(CLCCS(A, B, X, U), X0), Options(:T=>Inf))
 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)