Update transformation to canonical form

REQUIRE

@@ -3,7 +3,7 @@ Compat
 Expokit
 HybridSystems
 LazySets
-MathematicalSystems
+MathematicalSystems 0.6.3
 Memento
 Optim
 Plots

docs/src/lib/systems.md

@@ -48,3 +48,9 @@ vary over time.
 
 Time-varying nondeterministic inputs are chosen from a set of values that
 changes over time (with each time step).
+
+## Normalization
+
+```@docs
+normalize
+```

src/Utils/normalization.jl

@@ -1,67 +1,119 @@
-"""
-    normalize(system)
+# linear ivp in canonical form x' = Ax
+const LCF = LCS{N, AN} where {N, AN<:AbstractMatrix{N}}
 
-Normalize a continuous system.
+# 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}
 
-### Input
+# type union of canonical forms
+const CF = Union{<:LCF, <:ACF}
 
-- `system` -- continuous system
+# accepted types of non-deterministic inputs (non-canonical form)
+const UNCF = Union{<:LazySet{N}, Vector{<:LazySet{N}}, <:AbstractInput} where {N}
 
-### Output
+# accepted types for the state constraints (non-canonical form)
+const XNCF = Union{<:LazySet{N}, Nothing} where {N}
 
-Either the same system if it already conforms to our required structure, or a
-new system otherwise.
+"""
+    normalize(system::AbstractSystem)
 
-### Algorithm
+Transform a mathematical system to a normalized (or canonical) form.
 
-We apply the following normalization steps.
+### Input
 
-* [`normalize_inputs`](@ref)
-"""
-function normalize(system::InitialValueProblem{<:Union{AbstractContinuousSystem,
-                                                       AbstractDiscreteSystem}})
-    return normalize_inputs(system)
-end
+- `system` -- system; it can be discrete or continuous
 
-"""
-    normalize_inputs(system)
+### Output
 
-Normalize the inputs of a continuous system.
+Either the same system if it already conforms to a canonical form, or a new
+system otherwise.
 
-### Input
+### Notes
 
-- `system` -- continuous system
+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.
 
-### Output
+### Algorithm
 
-Either the same system if the inputs are of type `AbstractInput`, or a new
-system that wraps the inputs in a `ConstantInput`.
+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.
+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.
+
+The generalization to canonical systems with constraints and possibly time-varying
+non-deterministic inputs is considered next. These systems are of the form
+``x' = Ax + u, u ∈ \\mathcal{U}, x ∈ \\mathcal{X}``. The system type is
+`ConstrainedLinearControlContinuousSystem`, where `A` is a matrix, `X` is a set
+and `U` is an input, that is, any concrete subtype of `AbstractInput`.
+
+If `U` is not given as an input, normalization accepts either a `LazySet`, or
+a vector of `LazySet`s. In these cases, the sets are wrapped around an appropriate
+concrete input type.
+
+If the system does not conform to a canonical form, the implementation tries
+to make the transformation; otherwise an error is thrown. In particular, ODEs
+of the form ``x' = Ax + Bu`` are mapped into ``x' = Ax + u, u ∈ B\\mathcal{U}``,
+where now ``u`` has the same dimensions as ``x``.
+
+The transformations described above are analogous in the discrete case, i.e.
+``x_{k+1} = A x_k`` and ``x_{k+1} = Ax_{k} + u_k, u_k ∈ \\mathcal{U}, x_k ∈ \\mathcal{X}``
+for the linear and affine cases respectively.
 """
-function normalize_inputs(system)
-    inputdim(system) == 0 && return system
-    U = inputset(system)
-    U isa AbstractInput && return system
-    if !(U isa LazySet)
-        throw(ArgumentError("inputs of type $(typeof(U)) are not supported"))
-    end
-    return IVP(wrap_inputs(system.s, U), system.x0)
+function normalize(system::AbstractSystem)
+    throw(ArgumentError("the system type $(typeof(system)) is currently not supported"))
 end
 
-function wrap_inputs(system::CLCCS, U::LazySet)
-    return CLCCS(system.A, system.B, system.X, ConstantInput(U))
-end
-
-function wrap_inputs(system::CLCDS, U::LazySet)
-    return CLCDS(system.A, system.B, system.X, ConstantInput(U))
+# systems already in canonical form, i.e. which don't need normalization
+normalize(system::CF) = system
+
+# 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)
+    @eval begin
+        function normalize(system::$CLCS{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)
+            U = _wrap_inputs(system.U, system.B)
+            $CLCS(system.A, I(n, N), X, U)
+        end
+    end
 end
 
-function wrap_inputs(system::CACCS, U::LazySet)
-    return CACCS(system.A, system.B, system.c, system.X, ConstantInput(U))
+# 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)
+    @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}
+            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)
+        end
+    end
 end
 
-function wrap_inputs(system::CACDS, U::LazySet)
-    return CACDS(system.A, system.B, system.c, system.X, ConstantInput(U))
-end
+_wrap_invariant(X::LazySet, n::Int) = X
+_wrap_invariant(X::Nothing, n::Int) = Universe(n)
+
+_wrap_inputs(U::AbstractInput, B::IdentityMultiple) = U
+_wrap_inputs(U::AbstractInput, B::AbstractMatrix) = map(u -> B*u, U)
+_wrap_inputs(U::LazySet, B::IdentityMultiple) = ConstantInput(U)
+_wrap_inputs(U::LazySet, B::AbstractMatrix) = ConstantInput(B*U)
+_wrap_inputs(U::Vector{<:LazySet}, B::IdentityMultiple) = VaryingInput(U)
+_wrap_inputs(U::Vector{<:LazySet}, B::AbstractMatrix) = VaryingInput(map(u -> B*u, U))
+
+_wrap_inputs(U::AbstractInput, B::IdentityMultiple, c::AbstractVector) = U ⊕ c
+_wrap_inputs(U::AbstractInput, B::AbstractMatrix, c::AbstractVector) = map(u -> B*u ⊕ c, U)
+_wrap_inputs(U::LazySet, B::IdentityMultiple, c::AbstractVector) = ConstantInput(U ⊕ c)
+_wrap_inputs(U::LazySet, B::AbstractMatrix, c::AbstractVector) = ConstantInput(B*U ⊕ c)
+_wrap_inputs(U::Vector{<:LazySet}, B::IdentityMultiple, c::AbstractVector) = VaryingInput(map(u -> u ⊕ c, U))
+_wrap_inputs(U::Vector{<:LazySet}, B::AbstractMatrix, c::AbstractVector) = VaryingInput(map(u -> B*u ⊕ c, U))
 
 """
     distribute_initial_set(system::InitialValueProblem{<:HybridSystem, <:LazySet)

src/solve.jl

@@ -56,25 +56,26 @@ end
 solve(system::AbstractSystem, options::Pair{Symbol,<:Any}...) =
     solve(system, Options(Dict{Symbol,Any}(options)))
 
-function solve!(system::InitialValueProblem{<:Union{AbstractContinuousSystem,
-                                                     AbstractDiscreteSystem}},
+function solve!(problem::InitialValueProblem,
                 options_input::Options;
                 op::ContinuousPost=default_operator(system),
                 invariant::Union{LazySet, Nothing}=nothing
                )::AbstractSolution
-    system = normalize(system)
-    options = init!(op, system, options_input)
+
+    system = normalize(problem.s)
+    problem = IVP(system, problem.x0)
+    options = init!(op, problem, options_input)
 
     # coordinate transformation
     if options[:coordinate_transformation] != ""
         info("Transformation...")
-        (system, transformation_matrix) =
-            @timing transform(system, options[:coordinate_transformation])
+        (problem, transformation_matrix) =
+            @timing transform(problem, options[:coordinate_transformation])
         options[:transformation_matrix] = transformation_matrix
         invariant = options[:coordinate_transformation] * invariant
     end
 
-    post(op, system, invariant, options)
+    post(op, problem, invariant, options)
 end
 
 """

test/Reachability/solve_continuous.jl

@@ -120,7 +120,7 @@ s = solve(system, Options(:T=>0.1),
 # =======================================================
 # Affine ODE with a nondeterministic input, x' = Ax + Bu
 # =======================================================
-# linear ODE: 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
@@ -137,16 +137,13 @@ U = Interval(0.99, 1.01) × Interval(0.99, 1.01)
 # initial set
 X0 = BallInf(ones(4), 0.1)
 
-# dense identity matrix
-E = Matrix(1.0I, size(A))
-
 # inputs
-U1 = ConstantInput(B*U)
-U2 = B*U  # use internal conversion
+U1 = ConstantInput(U)
+U2 = U  # use internal wrapping
 
 for inputs in [U1, U2]
     # instantiate system
-    Δ = CLCCS(A, E, nothing, inputs)
+    Δ = CLCCS(A, B, nothing, U)
     𝒮 = InitialValueProblem(Δ, X0)
 
     # default options (computes all variables)