diff --git a/docs/src/lib/discretize.md b/docs/src/lib/discretize.md index a651ca3c..8f03ca7f 100644 --- a/docs/src/lib/discretize.md +++ b/docs/src/lib/discretize.md @@ -11,17 +11,7 @@ CurrentModule = Reachability.ReachSets ```@docs discretize -``` - -## Dense-time reachability - -```@docs -discr_bloat_firstorder -discr_bloat_interpolation -``` - -## Discrete-time reachability - -```@docs -discr_no_bloat +discretize_interpolation +discretize_firstorder +discretize_nobloating ``` diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index ac7a9560..94cc81ad 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -24,8 +24,7 @@ end function options_BFFPSV18() return OptionSpec[ # general options - OptionSpec(:approx_model, "forward", domain=String, domain_check=( - v -> v in ["forward", "backward", "firstorder", "nobloating"]), + OptionSpec(:discretization, "forward", domain=String, info="model for bloating/continuous time analysis"), OptionSpec(:algorithm, "explicit", domain=String, domain_check=( v -> v in ["explicit", "wrap"]), info="algorithm backend"), @@ -41,16 +40,10 @@ function options_BFFPSV18() "containing its indices"), # discretization options - OptionSpec(:lazy_sih, false, domain=Bool, - info="use a lazy symmetric interval hull in discretization?"), - OptionSpec(:lazy_expm, false, domain=Bool, - info="use a lazy matrix exponential all the time?"), - OptionSpec(:lazy_expm_discretize, false, domain=Bool, - info="use a lazy matrix exponential in discretization?"), - OptionSpec(:pade_expm, false, domain=Bool, - info="use the Padé approximant method (instead of Julia's " * - " built-in 'exp') to compute the lazy matrix exponential " * - "in discretization?"), + OptionSpec(:sih_method, "concrete", domain=String, + info="method to compute the symmetric interval hull in discretization"), + OptionSpec(:exp_method, "base", domain=String, + info="method to compute the matrix exponential"), OptionSpec(:assume_sparse, false, domain=Bool, info="use an analysis for sparse discretized matrices?"), @@ -208,12 +201,6 @@ function normalization_BFFPSV18!(𝑂::TwoLayerOptions) end function validation_BFFPSV18(𝑂) - # lazy_expm_discretize & lazy_expm - if !𝑂[:lazy_expm_discretize] && 𝑂[:lazy_expm] - throw(DomainError(𝑂[:lazy_expm_discretize], "cannot use option " * - "':lazy_expm' with deactivated option ':lazy_expm_discretize'")) - end - # block_types if haskey_specified(𝑂, :block_types) for (key, value) in 𝑂[:block_types] @@ -337,7 +324,6 @@ Calculate the reachable states of the given initial value problem using `BFFPSV1 - `𝑂` -- algorithm-specific options """ function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, invariant, 𝑂::Options) - # TODO temporary hack for refactoring 𝑂 = TwoLayerOptions(merge(𝑂, 𝒫.options.specified), 𝒫.options.defaults) # convert matrix diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index b0356d2b..2736a12b 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -182,16 +182,9 @@ function check_property(S::IVP{<:AbstractContinuousSystem}, # Time discretization # =================== info("Time discretization...") - Δ = @timing begin - discretize( - S, - options[:δ], - approx_model=options[:approx_model], - pade_expm=options[:pade_expm], - lazy_expm=options[:lazy_expm_discretize], - lazy_sih=options[:lazy_sih] - ) - end - Δ = matrix_conversion_lazy_explicit(Δ, options) + Δ = @timing discretize(S, options[:δ], algorithm=options[:discretization], + exp_method=options[:exp_method], + sih_method=options[:sih_method]) + Δ = matrix_conversion(Δ, options) return check_property(Δ, property, options) end diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl index 3ce6d5eb..9ee0e64c 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -227,17 +227,11 @@ function reach(system::IVP{<:AbstractContinuousSystem}, # Time discretization # =================== info("Time discretization...") - Δ = @timing begin - discretize( - system, - options[:δ], - approx_model=options[:approx_model], - pade_expm=options[:pade_expm], - lazy_expm=options[:lazy_expm_discretize], - lazy_sih=options[:lazy_sih] - ) - end - Δ = matrix_conversion_lazy_explicit(Δ, options) + Δ = @timing discretize(system, options[:δ], algorithm=options[:discretization], + exp_method=options[:exp_method], + sih_method=options[:sih_method]) + + Δ = matrix_conversion(Δ, options) return reach(Δ, invariant, options) end diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 23223202..dbdc0256 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -1,310 +1,652 @@ +const LDS = LinearDiscreteSystem +const CLCDS = ConstrainedLinearControlDiscreteSystem + +@inline Id(n) = Matrix(1.0I, n, n) + """ - discretize(cont_sys, δ; [approx_model], [pade_expm], [lazy_expm], [lazy_sih]) + discretize(𝑆, δ; [algorithm], [exp_method], [sih_method]) + +Apply an approximation model to `S` obtaining a discrete initial value problem. + +### Input + +- `𝑆` -- initial value problem for a continuous affine ODE with + non-deterministic inputs +- `δ` -- step size +- `algorithm` -- (optional, default: `"forward"`) the algorithm used to + compute the approximation model for the discretization, + choose among: + + - `"forward"` -- use forward-time interpolation + - `"backward"` -- use backward-time interpolation + - `"firstorder"` -- use first-order approximation of the ODE + - `"nobloating"` -- do not bloat the initial states + +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: + + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` + +- `sih_method` -- (optional, default: `"concrete"`) the method used to take the + symmetric interval hull operation, choose among: + + - `"concrete"` -- compute the full symmetric interval hull using the function + `symmetric_interval_hull` from `LazySets.Approximations` + - `"lazy"` -- compute a wrapper set type around symmetric interval hull + in a lazy way using `SymmetricIntervalHull` -Discretize a continuous system of ODEs with nondeterministic inputs. +### Output -## Input +The initial value problem of a discrete system. -- `cont_sys` -- continuous system -- `δ` -- step size -- `approx_model` -- the method to compute the approximation model for the - discretization, among: +### Algorithm - - `forward` -- use forward-time interpolation - - `backward` -- use backward-time interpolation - - `firstorder` -- use first order approximation of the ODE - - `nobloating` -- do not bloat the initial states - (use for discrete-time reachability) +Let ``𝑆 : x' = Ax(t) + u(t)``, ``x(0) ∈ \\mathcal{X}_0``, ``u(t) ∈ U`` be the +given continuous affine ODE `𝑆`, where `U` is the set of non-deterministic inputs +and ``\\mathcal{X}_0`` is the set of initial states. Recall that the system +`𝑆` is called homogeneous whenever `U` is the empty set. -- `pade_expm` -- (optional, default = `false`) if true, use Pade approximant - method to compute matrix exponentials of sparse matrices; - otherwise use Julia's buil-in `expm` -- `lazy_expm` -- (optional, default = `false`) if true, compute the matrix - exponential in a lazy way (suitable for very large systems) -- `lazy_sih` -- (optional, default = `true`) if true, compute the - symmetric interval hull in a lazy way (suitable if only a - few dimensions are of interest) +Given a step size ``δ``, this function computes a set, `Ω₀`, that guarantees to +contain all the trajectories of ``𝑆`` starting at any ``x(0) ∈ \\mathcal{X}_0`` +and for any input function that satisfies ``u(t) ∈ U``, for any ``t ∈ [0, δ]``. -## Output +The initial value problem returned by this function consists of the set `Ω₀` +together with the coefficient matrix ``ϕ = e^{Aδ}`` and a transformed +set of inputs if `U` is non-empty. -A discrete system. +In the literature, the method to obtain `Ω₀` is called the *approximation model* +and different alternatives have been proposed. See the argument `algorithm` +for available options. For the reference to the original papers, see the +docstring of each method `discretize_...`. -## Notes +In the dense-time case, the transformation is such that the trajectories +of the given continuous system are included in the computed flowpipe of the +discretized system. -This function applies an approximation model to transform a continuous affine -system into a discrete affine system. -This transformation allows to do dense time reachability, i.e. such that the -trajectories of the given continuous system are included in the computed -flowpipe of the discretized system. -For discrete-time reachability, use `approx_model="nobloating"`. +In the discrete-time case, there is no bloating of the initial states and the +input is assumed to remain constant between sampled times. Use the option +`algorithm="nobloating"` for this setting. + +Several methods to compute the matrix exponential are availabe. Use `exp_method` +to select one. For very large systems, computing the full matrix exponential is +expensive hence it is preferable to compute the action of the matrix exponential +over vectors when needed, `e^{δA} v` for each `v`. Use the option +`exp_method="lazy"` for this purpose. """ -function discretize(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, +function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; - approx_model::String="forward", - pade_expm::Bool=false, - lazy_expm::Bool=false, - lazy_sih::Bool=true)::InitialValueProblem{<:AbstractDiscreteSystem} - - if approx_model in ["forward", "backward"] - return discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm, - lazy_expm, lazy_sih) - elseif approx_model == "firstorder" - return discr_bloat_firstorder(cont_sys, δ) - elseif approx_model == "nobloating" - return discr_no_bloat(cont_sys, δ, pade_expm, lazy_expm) + algorithm::String="forward", + exp_method::String="base", + sih_method::String="concrete") + + if algorithm in ["forward", "backward"] + return discretize_interpolation(𝑆, δ, algorithm=algorithm, + exp_method=exp_method, sih_method=sih_method) + elseif algorithm == "firstorder" + return discretize_firstorder(𝑆, δ, exp_method=exp_method) + elseif algorithm == "nobloating" + return discretize_nobloating(𝑆, δ, exp_method=exp_method) else - error("The approximation model is invalid") + throw(ArgumentError("the algorithm $algorithm is unknown")) end end """ - bloat_firstorder(cont_sys, δ) + exp_Aδ(A::AbstractMatrix, δ::Float64; [exp_method]) -Compute bloating factors using first order approximation. +Compute the matrix exponential ``e^{Aδ}``. -## Input +### Input -- `cont_sys` -- a continuous affine system -- `δ` -- step size +- `A` -- coefficient matrix +- `δ` -- step size +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: -## Notes + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` -In this algorithm, the infinity norm is used. -See also: `discr_bloat_interpolation` for more accurate (less conservative) -bounds. +### Output -## Algorithm +A matrix. +""" +function exp_Aδ(A::AbstractMatrix{Float64}, δ::Float64; exp_method="base") + if exp_method == "base" + return expmat(Matrix(A*δ)) + elseif exp_method == "lazy" + return SparseMatrixExp(A*δ) + elseif exp_method == "pade" + return padm(A*δ) + else + throw(ArgumentError("the exponentiation method $exp_method is unknown")) + end +end -This uses a first order approximation of the ODE, and matrix norm upper bounds, -see Le Guernic, C., & Girard, A., 2010, *Reachability analysis of linear systems -using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.* +@inline function Pmatrix(A, δ, n) + return [A*δ sparse(δ*I, n, n) spzeros(n, n) ; + spzeros(n, 2*n ) sparse(δ*I, n, n); + spzeros(n, 3*n ) ] +end + +""" + Φ₁(A, δ; [exp_method]) + +Compute the series + +```math +Φ₁(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+1}}{(i+1)!}A^i, +``` +where ``A`` is a square matrix of order ``n`` and ``δ ∈ \\mathbb{R}_{≥0}``. + +### Input + +- `A` -- coefficient matrix +- `δ` -- step size +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: + + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` +### Output + +A matrix. + +### Algorithm + +We use the method from [1]. If ``A`` is invertible, ``Φ₁`` can be computed as + +```math +Φ₁(A, δ) = A^{-1}(e^{δA} - I_n). +``` + +In the general case, implemented in this function, it can be computed as +submatrices of the block matrix + +```math +P = \\exp \\begin{pmatrix} +Aδ && δI_n && 0 \\ +0 && 0 && δI_n \\ +0 && 0 && 0 +\\end{array}. +``` +It can be shown that `Φ₁(A, δ) = P[1:n, (n+1):2*n]`. + +[1] Frehse, Goran, et al. "SpaceEx: Scalable verification of hybrid systems." +International Conference on Computer Aided Verification. Springer, Berlin, +Heidelberg, 2011. """ -function discr_bloat_firstorder(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, - δ::Float64) - - A, X0 = cont_sys.s.A, cont_sys.x0 - Anorm = norm(Matrix(A), Inf) - ϕ = expmat(Matrix(A)) - RX0 = norm(X0, Inf) - - if inputdim(cont_sys) == 0 - # linear case - α = (exp(δ*Anorm) - 1. - δ*Anorm) * RX0 - Ω0 = CH(X0, ϕ * X0 + Ball2(zeros(size(ϕ, 1)), α)) - return DiscreteSystem(ϕ, Ω0) +function Φ₁(A, δ; exp_method="base") + n = size(A, 1) + if exp_method == "base" + P = expmat(Matrix(Pmatrix(A, δ, n))) + Φ₁_Aδ = P[1:n, (n+1):2*n] + + elseif exp_method == "lazy" + P = SparseMatrixExp(Pmatrix(A, δ, n)) + Φ₁_Aδ = sparse(get_columns(P, (n+1):2*n)[1:n, :]) + + elseif exp_method == "pade" + P = padm(Pmatrix(A, δ, n)) + Φ₁_Aδ = P[1:n, (n+1):2*n] + else - # affine case; TODO: unify Constant and Varying input branches? - Uset = inputset(cont_sys) + throw(ArgumentError("the exponentiation method $exp_method is unknown")) + end + + return Φ₁_Aδ +end + +""" + Φ₂(A, δ; [exp_method]) + +Compute the series + +```math +Φ₂(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+2}}{(i+2)!}A^i, +``` +where ``A`` is a square matrix of order ``n`` and ``δ ∈ \\mathbb{R}_{≥0}``. + +### Input + +- `A` -- coefficient matrix +- `δ` -- step size +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: + + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` + +### Output + +A matrix. + +### Algorithm + +We use the method from [1]. If ``A`` is invertible, ``Φ₂`` can be computed as + +```math +Φ₂(A, δ) = A^{-2}(e^{δA} - I_n - δA). +``` + +In the general case, implemented in this function, it can be computed as +submatrices of the block matrix + +```math +P = \\exp \\begin{pmatrix} +Aδ && δI_n && 0 \\ +0 && 0 && δI_n \\ +0 && 0 && 0 +\\end{array}. +``` +It can be shown that `Φ₂(A, δ) = P[1:n, (2*n+1):3*n]`. + +[1] Frehse, Goran, et al. "SpaceEx: Scalable verification of hybrid systems." +International Conference on Computer Aided Verification. Springer, Berlin, +Heidelberg, 2011. +""" +function Φ₂(A, δ; exp_method="base") + n = size(A, 1) + if exp_method == "base" + P = expmat(Matrix(Pmatrix(A, δ, n))) + Φ₂_Aδ = P[1:n, (2*n+1):3*n] + + elseif exp_method == "lazy" + P = SparseMatrixExp(Pmatrix(A, δ, n)) + Φ₂_Aδ = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) + + elseif exp_method == "pade" + P = padm(Pmatrix(A, δ, n)) + Φ₂_Aδ = P[1:n, (2*n+1):3*n] + + else + throw(ArgumentError("the exponentiation method $exp_method is unknown")) + end + + return Φ₂_Aδ +end + +""" + discretize_firstorder(𝑆, δ; [p], [exp_method]) + +Apply a first-order approximation model to `S` obtaining a discrete initial +value problem. + +### Input + +- `𝑆` -- initial value problem for a continuous affine ODE with + non-deterministic inputs +- `δ` -- step size +- `p` -- (optional, default: `Inf`) parameter in the considered norm +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: + + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` + +### Output + +The initial value problem for a discrete system. In particular: + +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. + +### Algorithm + +Let us define some notation. Let + +```math +𝑆 : x' = Ax(t) + u(t) +``` +with ``x(0) ∈ \\mathcal{X}_0``, ``u(t) ∈ U`` be the given continuous affine ODE +`𝑆`, where `U` is the set of non-deterministic inputs and ``\\mathcal{X}_0`` +is the set of initial states. + +Define ``R_{\\mathcal{X}_0} = \\max_{x ∈ \\mathcal{X}_0} ‖x‖``, +`D_{\\mathcal{X}_0} = \\max_{x, y ∈ \\mathcal{X}_0} ‖x-y‖`` and +``R_{U} = \\max_{u ∈ U} ‖u‖``. If only the support functions of ``\\mathcal{X}_0`` +and ``U`` are known, these values might be hard to compute for some norms. See +`Notes` below for details. + +Let ``Ω₀`` be the set defined as: +```math +Ω₀ = ConvexHull(\\mathcal{X}_0, e^{δA}\\mathcal{X}_0 ⊕ δU ⊕ αB_p) +``` +where ``α = (e^{δ ‖A‖} - 1 - δ‖A‖)*(R_{\\mathcal{X}_0} + R_{U} / ‖A‖)`` and +``B_p`` denotes the unit ball for the considered ``p``-norm. + +It is proved in [Lemma 1, 1] that the set of states reachable by ``S`` in the time +interval ``[0, δ]``, which we denote ``R_{[0,δ]}(\\mathcal{X}_0)`` here, +is included in ``Ω₀``: + +```math +R_{[0,δ]}(\\mathcal{X}_0) ⊆ Ω₀. +``` + +Moreover, if ``d_H(A, B)`` denotes the Hausdorff distance between the sets ``A`` +and ``B`` in ``\\mathbb{R}^n``, then + +```math +d_H(Ω₀, R_{[0,δ]}(\\mathcal{X}_0)) ≤ \\frac{1}{4}(e^{δ ‖A‖} - 1) D_{\\mathcal{X}_0} + 2α. +``` +Hence, the approximation error can be made arbitrarily small by choosing ``δ`` +small enough. + +Here we allow ``U`` to be a sequence of time varying non-deterministic inputs. + +### Notes + +In this implementation, the infinity norm is used by default. Other usual norms +are ``p=1`` and ``p=2``. However, note that not all norms are supported; see the +documentation of `?norm` in `LazySets` for the supported norms. + +See also [`discretize_interpolation`](@ref) for an alternative algorithm that +uses less conservative bounds. + +[1] Le Guernic, C., & Girard, A., 2010, *Reachability analysis of linear systems +using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.* +""" +function discretize_firstorder(𝑆::InitialValueProblem, + δ::Float64; + p::Float64=Inf, + exp_method::String="base") + + # unwrap coefficient matrix and initial states + A, X0 = 𝑆.s.A, 𝑆.x0 + + # system size; A is assumed square + n = size(A, 1) + + Anorm = norm(Matrix(A), p) + κ = exp(δ*Anorm) - 1.0 - δ*Anorm + RX0 = norm(X0, p) + + # compute exp(A*δ) + ϕ = exp_Aδ(A, δ, exp_method=exp_method) + + if inputdim(𝑆) == 0 + α = κ * RX0 + □α = Ballp(p, zeros(n), α) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ □α) + return IVP(LDS(ϕ), Ω0) + + elseif isaffine(𝑆) + Uset = inputset(𝑆) if Uset isa ConstantInput U = next_set(Uset) - RU = norm(U, Inf) - α = (exp(δ*Anorm) - 1. - δ*Anorm)*(RX0 + RU/Anorm) - β = (exp(δ*Anorm) - 1. - δ*Anorm)*RU/Anorm - Ω0 = CH(X0, ϕ * X0 + δ * U + Ball2(zeros(size(ϕ, 1)), α)) - discr_U = δ * U + Ball2(zeros(size(ϕ, 1)), β) - return DiscreteSystem(ϕ, Ω0, discr_U) + RU = norm(U, p) + + # bloating coefficients + α = κ*(RX0 + RU/Anorm) + β = κ*RU/Anorm + + # transformation of the initial states + □α = Ballp(p, zeros(n), α) + Ω0 = ConvexHull(X0, ϕ*X0 ⊕ δ*U ⊕ □α) + + # transformation of the inputs + □β = Ballp(p, zeros(n), β) + Ud = ConstantInput(δ*U ⊕ □β) + return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) + elseif Uset isa VaryingInput - discr_U = Vector{LazySet}(undef, length(Uset)) + Ud = Vector{LazySet}(undef, length(Uset)) for (i, Ui) in enumerate(Uset) - RU = norm(Ui, Inf) - α = (exp(δ*Anorm) - 1. - δ*Anorm)*(RX0 + RU/Anorm) - β = (exp(δ*Anorm) - 1. - δ*Anorm)*RU/Anorm - Ω0 = CH(X0, ϕ * X0 + δ * Ui + Ball2(zeros(size(ϕ, 1)), α)) - discr_U[i] = δ * Ui + Ball2(zeros(size(ϕ, 1)), β) + RU = norm(Ui, p) + + # bloating factors + α = κ*(RX0 + RU/Anorm) + β = κ*RU/Anorm + + if i == 1 + # transform initial states + □α = Ballp(p, zeros(n), α) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ * Ui ⊕ □α) + end + + # transform inputs + □β = Ballp(p, zeros(n), β) + Ud[i] = δ*Ui ⊕ □β end - return DiscreteSystem(ϕ, Ω0, discr_U) + Ud = VaryingInput(Ud) + return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) + else + throw(ArgumentError("input of type $(typeof(U)) is not allowed")) end + else + throw(ArgumentError("this function only applies to linear or affine systems")) end - - end """ - discr_no_bloat(cont_sys, δ, pade_expm, lazy_expm) + discretize_nobloating(𝑆, δ; [exp_method]) Discretize a continuous system without bloating of the initial states, suitable for discrete-time reachability. -## Input +### Input -- `cont_sys` -- a continuous system -- `δ` -- step size -- `pade_expm` -- if `true`, use Pade approximant method to compute the - matrix exponential -- `lazy_expm` -- if `true`, compute the matrix exponential in a lazy way - (suitable for very large systems) +- `𝑆` -- a continuous system +- `δ` -- step size +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: -## Output + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` -A discrete system. +### Output -## Algorithm +The initial value problem for a discrete system. In particular: + +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. + +### Algorithm -The transformation implemented here is the following: +Let us define some notation. Let -- `A -> Phi := exp(A*delta)` -- `U -> V := M*U` -- `X0 -> X0hat := X0` +```math +𝑆 : x' = Ax(t) + u(t) +``` +with ``x(0) ∈ \\mathcal{X}_0``, ``u(t) ∈ U`` be the given continuous affine ODE +`𝑆`, where `U` is the set of non-deterministic inputs and ``\\mathcal{X}_0`` +is the set of initial states. -where `M` corresponds to `Phi1(A, delta)` in Eq. (8) of *SpaceEx: Scalable -Verification of Hybrid Systems.* +The approximation model implemented in this function is such that there is no bloating, +i.e. we don't bloat the initial states and don't multiply the input by the step +size δ, as required for the dense time case. -In particular, there is no bloating, i.e. we don't bloat the initial states and -dont multiply the input by the step size δ, as required for the dense time case. +The transformations are: + +- ``Φ ← \\exp^{Aδ}`` +- ``Ω₀ ← \\mathcal{X}_0`` +- ``V ← Φ₁(A, δ)U(k)``, where ``Φ₁(A, δ)`` is defined in + [`Φ₁(A, δ; [exp_method])`](@ref). + +Here we allow ``U`` to be a sequence of time varying non-deterministic input sets. """ -function discr_no_bloat(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, - δ::Float64, - pade_expm::Bool, - lazy_expm::Bool) +function discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, + δ::Float64; + exp_method::String="base") - A, X0 = cont_sys.s.A, cont_sys.x0 - n = size(A, 1) - if lazy_expm - ϕ = SparseMatrixExp(A * δ) - else - if pade_expm - ϕ = padm(A * δ) - else - ϕ = expmat(Matrix(A * δ)) - end - end + # unwrap coefficient matrix and initial states + A, X0 = 𝑆.s.A, 𝑆.x0 + + # compute matrix ϕ = exp(Aδ) + ϕ = exp_Aδ(A, δ, exp_method=exp_method) + + # initial states remain unchanged + Ω0 = copy(X0) # early return for homogeneous systems - if cont_sys isa IVP{<:LinearContinuousSystem} - Ω0 = X0 - return DiscreteSystem(ϕ, Ω0) + if inputdim(𝑆) == 0 + return IVP(LDS(ϕ), Ω0) end - U = inputset(cont_sys) - inputs = next_set(U, 1) # compute matrix to transform the inputs - if lazy_expm - P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) - Phi1Adelta = sparse(get_columns(P, (n+1):2*n)[1:n, :]) - else - if pade_expm - P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) - else - P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)])) - end - Phi1Adelta = P[1:n, (n+1):2*n] - end - - discretized_U = Phi1Adelta * inputs + Phi1Adelta = Φ₁(A, δ, exp_method=exp_method) + U = inputset(𝑆) + Ud = map(ui -> Phi1Adelta * ui, U) - Ω0 = X0 - - if U isa ConstantInput - return DiscreteSystem(ϕ, Ω0, discretized_U) - else - discretized_U = VaryingInput([Phi1Adelta * Ui for Ui in U]) - return DiscreteSystem(ϕ, Ω0, discretized_U) - end + return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0) end """ - discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm, lazy_expm) + discretize_interpolation(𝑆, δ; [algorithm], [exp_method], [sih_method]) Compute bloating factors using forward or backward interpolation. -## Input +### Input + +- `𝑆` -- a continuous system +- `δ` -- step size +- `algorithm` -- choose the algorithm to compute the approximation model + among `"forward"` and `"backward"` +- `exp_method` -- (optional, default: `"base"`) the method used to take the matrix + exponential of the coefficient matrix, choose among: -- `cs` -- a continuous system -- `δ` -- step size -- `approx_model` -- choose the approximation model among `"forward"` and - `"backward"` -- `pade_expm` -- if true, use Pade approximant method to compute the - matrix exponential -- `lazy_expm` -- if true, compute the matrix exponential in a lazy way - suitable for very large systems) + - `"base"` -- the scaling and squaring method implemented in Julia base, + see `?exp` for details + - `"pade"` -- use Pade approximant method to compute matrix exponentials of + sparse matrices, implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` + +- `sih_method` -- (optional, default: `"concrete"`) the method used to take the + symmetric interval hull operation, choose among: + + - `"concrete"` -- compute the full symmetric interval hull + - `"lazy"` -- compute a wrapper set type around symmetric interval hull in a + lazy way + +### Output + +The initial value problem for a discrete system. In particular: + +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ## Algorithm -See Frehse et al., CAV'11, *SpaceEx: Scalable Verification of Hybrid Systems*, -Lemma 3. +Let us define some notation. Let -Note that in the unlikely case that A is invertible, the result can also -be obtained directly, as a function of the inverse of A and `e^{At} - I`. +```math +𝑆 : x' = Ax(t) + u(t) +``` +with ``x(0) ∈ \\mathcal{X}_0``, ``u(t) ∈ U`` be the given continuous affine ODE +`𝑆`, where `U` is the set of non-deterministic inputs and ``\\mathcal{X}_0`` +is the set of initial states. -The matrix `P` is such that: `ϕAabs = P[1:n, 1:n]`, -`Phi1Aabsdelta = P[1:n, (n+1):2*n]`, and `Phi2Aabs = P[1:n, (2*n+1):3*n]`. -""" -function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, - δ::Float64, - approx_model::String, - pade_expm::Bool, - lazy_expm::Bool, - lazy_sih::Bool) +The transformations are: - sih = lazy_sih ? SymmetricIntervalHull : symmetric_interval_hull +- ``Φ ← \\exp^{Aδ}``, +- ``Ω₀ ← ConvexHull(\\mathcal{X}_0, Φ\\mathcal{X}_0 ⊕ δU(0) ⊕ Eψ(U(0), δ) ⊕ E^+(\\mathcal{X}_0, δ))``, +- ``V ← δU(k) ⊕ Eψ(U(k), δ)``. - A, X0 = cont_sys.s.A, cont_sys.x0 - n = size(A, 1) +Here we allow ``U`` to be a sequence of time varying non-deterministic input sets. - # compute matrix ϕ = exp(Aδ) - if lazy_expm - ϕ = SparseMatrixExp(A*δ) +For the definition of the sets ``Eψ`` and ``E^+`` see [1]. The "backward" method +uses ``E^-``. + +[1] Frehse, Goran, et al. "SpaceEx: Scalable verification of hybrid systems." +International Conference on Computer Aided Verification. Springer, Berlin, +Heidelberg, 2011. +""" +function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, + δ::Float64; + algorithm::String="forward", + exp_method::String="base", + sih_method::String="concrete") + + if sih_method == "concrete" + sih = symmetric_interval_hull + elseif sih_method == "lazy" + sih = SymmetricIntervalHull else - if pade_expm - ϕ = padm(A*δ) - else - ϕ = expmat(Matrix(A*δ)) - end + throw(ArgumentError("the method $sih_method is unknown")) end - # early return for homogeneous systems - if cont_sys isa IVP{<:LinearContinuousSystem} - Ω0 = CH(X0, ϕ * X0) - return DiscreteSystem(ϕ, Ω0) - end - U = inputset(cont_sys) - inputs = next_set(U, 1) + # unwrap coefficient matrix and initial states + A, X0 = 𝑆.s.A, 𝑆.x0 + + # compute matrix ϕ = exp(Aδ) + ϕ = exp_Aδ(A, δ, exp_method=exp_method) # compute the transformation matrix to bloat the initial states - if lazy_expm - P = SparseMatrixExp([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) - Phi2Aabs = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) + Phi2Aabs = Φ₂(abs.(A), δ, exp_method=exp_method) + + if algorithm == "forward" + Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use E⁺ + elseif algorithm == "backward" + Einit = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) # use E⁻ else - if pade_expm - P = padm([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) - else - P = expmat(Matrix([abs.(A*δ) sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)])) - end - Phi2Aabs = P[1:n, (2*n+1):3*n] + throw(ArgumentError("the algorithm $approximation is unknown")) end - if isa(inputs, ZeroSet) - if approx_model == "forward" || approx_model == "backward" - Ω0 = CH(X0, ϕ * X0 + δ * inputs) - end - else - EPsi = sih(Phi2Aabs * sih(A * inputs)) - discretized_U = δ * inputs + EPsi - if approx_model == "forward" - EOmegaPlus = sih(Phi2Aabs * sih((A * A) * X0)) - Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaPlus) - elseif approx_model == "backward" - EOmegaMinus = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) - Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaMinus) - end + # early return for homogeneous systems + if inputdim(𝑆) == 0 + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ Einit) + return IVP(LDS(ϕ), Ω0) end + U = inputset(𝑆) + U0 = next_set(U, 1) + + Eψ0 = sih(Phi2Aabs * sih(A * U0)) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ*U0 ⊕ Eψ0 ⊕ Einit) + if U isa ConstantInput - return DiscreteSystem(ϕ, Ω0, discretized_U) + Ud = ConstantInput(δ*U0 ⊕ Eψ0) + + elseif U isa VaryingInput + Ud = Vector{LazySet}(undef, length(U)) + for (k, Uk) in enumerate(U) + Eψk = sih(Phi2Aabs * sih(A * Uk)) + Ud[k] = δ * Uk ⊕ Eψk + end + Ud = VaryingInput(Ud) + else - discretized_U = [δ * Ui + sih(Phi2Aabs * sih(A * Ui)) for Ui in U] - return DiscreteSystem(ϕ, Ω0, discretized_U) + throw(ArgumentError("input of type $(typeof(U)) is not allowed")) end + + return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0) end diff --git a/src/Reachability.jl b/src/Reachability.jl index 08ae9dd2..09efb06a 100644 --- a/src/Reachability.jl +++ b/src/Reachability.jl @@ -11,6 +11,12 @@ 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("Utils/Utils.jl") include("Options/dictionary.jl") diff --git a/src/Transformations/Transformations.jl b/src/Transformations/Transformations.jl index 5c7a16f3..789ffea4 100644 --- a/src/Transformations/Transformations.jl +++ b/src/Transformations/Transformations.jl @@ -84,10 +84,13 @@ function schur_transform(S::InitialValueProblem) T_inverse = F[:vectors] if S.s isa DiscreteSystem - return (DiscreteSystem(A_new, X0_new, U_new), T_inverse) + s = ConstrainedLinearControlDiscreteSystem(A_new, Matrix(1.0I, size(A_new)), nothing, U_new) + p = InitialValueProblem(s, X0_new) elseif S.s isa ContinuousSystem - return (ContinuousSystem(A_new, X0_new, U_new), T_inverse) + s = ConstrainedLinearControlContinuousSystem(A_new, Matrix(1.0I, size(A_new)), nothing, U_new) + p = InitialValueProblem(s, X0_new) end + return (p, T_inverse) end end # module diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 6d887a39..99da711e 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -508,34 +508,12 @@ function matrix_conversion(Δ, options; A_passed=nothing) if create_new_system # set new matrix if hasmethod(inputset, Tuple{typeof(Δ.s)}) - Δ = DiscreteSystem(A_new, Δ.x0, inputset(Δ)) + Δ = IVP(CLCDS(A_new, Matrix(1.0I, size(A_new)), nothing, inputset(Δ)), Δ.x0) else - Δ = DiscreteSystem(A_new, Δ.x0) + Δ = IVP(LDS(A_new), Δ.x0) end end return Δ end -# convert SparseMatrixExp to eplicit matrix -function matrix_conversion_lazy_explicit(Δ, options) - A = Δ.s.A - if !options[:lazy_expm] && options[:lazy_expm_discretize] - info("Making lazy matrix exponential explicit...") - @timing begin - n = options[:n] - if options[:assume_sparse] - B = sparse(Int[], Int[], eltype(A)[], n, n) - else - B = Matrix{eltype(A)}(n, n) - end - for i in 1:n - B[i, :] = get_row(A, i) - end - end - else - B = nothing - end - return matrix_conversion(Δ, options; A_passed=B) -end - end # module diff --git a/src/Utils/systems.jl b/src/Utils/systems.jl index 38cefd72..c5910317 100644 --- a/src/Utils/systems.jl +++ b/src/Utils/systems.jl @@ -15,6 +15,8 @@ import LazySets.constrained_dimensions *(M::AbstractMatrix, input::ConstantInput) = ConstantInput(M * input.U) +#= +# TODO: kept for backwards-compatibility. To be removed. # no input: x' = Ax, x(0) = X0 ContinuousSystem(A::AbstractMatrix, X0::LazySet) = IVP(LCS(A), X0) DiscreteSystem(A::AbstractMatrix, X0::LazySet) = IVP(LDS(A), X0) @@ -32,6 +34,7 @@ ContinuousSystem(A::AbstractMatrix, X0::LazySet, U::Vector{<:LazySet}) = Continu DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::VaryingInput) = IVP(CLCDS(A, convert(typeof(A), convert(typeof(A), Matrix{eltype(A)}(I, size(A)))), nothing, U), X0) DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::Vector{<:LazySet}) = DiscreteSystem(A, X0, VaryingInput(U)) +=# # convenience functions next_set(inputs::ConstantInput) = collect(nextinput(inputs, 1))[1] @@ -145,7 +148,7 @@ Adds an extra dimension to a continuous system. julia> using MathematicalSystems julia> A = sprandn(3, 3, 0.5); julia> X0 = BallInf(zeros(3), 1.0); -julia> s = ContinuousSystem(A, X0); +julia> s = InitialValueProblem(LinearContinuousSystem(A), X0); julia> sext = add_dimension(s); julia> statedim(sext) 4 @@ -154,8 +157,8 @@ julia> statedim(sext) If there is an input set, it is also extended: ```jldoctest add_dimension_cont_sys -julia> U = Ball2(ones(3), 0.1); -julia> s = ContinuousSystem(A, X0, U); +julia> U = ConstantInput(Ball2(ones(3), 0.1)); +julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0); julia> sext = add_dimension(s); julia> statedim(sext) 4 @@ -168,8 +171,8 @@ Extending a system with a varying input set: If there is an input set, it is also extended: ```jldoctest add_dimension_cont_sys -julia> U = [Ball2(ones(3), 0.1 * i) for i in 1:3]; -julia> s = ContinuousSystem(A, X0, U); +julia> U = VaryingInput([Ball2(ones(3), 0.1 * i) for i in 1:3]); +julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0); julia> sext = add_dimension(s); julia> statedim(sext) 4 @@ -192,10 +195,11 @@ function add_dimension(cs, m=1) X0ext = add_dimension(cs.x0, m) if hasmethod(inputset, Tuple{typeof(cs.s)}) Uext = map(x -> add_dimension(x, m), inputset(cs)) - return ContinuousSystem(Aext, X0ext, Uext) + s = ConstrainedLinearControlContinuousSystem(Aext, Matrix(1.0I, size(Aext)), nothing, Uext) else - return ContinuousSystem(Aext, X0ext) + s = LinearContinuousSystem(Aext) end + return InitialValueProblem(s, X0ext) end """ diff --git a/src/solve.jl b/src/solve.jl index d1730d6d..fcbd66ff 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -184,7 +184,7 @@ function solve!(system::InitialValueProblem{<:HybridSystem, # TODO temporary hack, to be resolved in #447 options_copy[:mode] = "reach" end - reach_tube = solve!(ContinuousSystem(loc.A, X0.X, loc.U), + reach_tube = solve!(IVP(loc, X0.X), options_copy, op=opC, invariant=source_invariant) diff --git a/test/Properties/unit_tune.jl b/test/Properties/unit_tune.jl index f53a3465..2b9507d9 100644 --- a/test/Properties/unit_tune.jl +++ b/test/Properties/unit_tune.jl @@ -5,8 +5,8 @@ We test the line plot!(x->x, x->-24*x+375, 0., 50., line=2., color="red", linest =# A = [0. 0.5 0. 0. ; 0. 0. 0. 0. ; 0. 0. 0. 0.7 ; 0. 0. 0. 0.] X0 = Singleton([0.,5.,100.,0]) -U = Singleton([0.,0.,0.,-9.81]) -S = ContinuousSystem(A, X0, U) +U = ConstantInput(Singleton([0.,0.,0.,-9.81])) +S = IVP(CLCDS(A, Matrix(1.0I, 4, 4), nothing, U), X0) time_horizon = 20.0 prec = 1e-4 initial_δ = 0.5 diff --git a/test/ReachSets/unit_discretization.jl b/test/ReachSets/unit_discretization.jl index 77297ea4..560bbae1 100644 --- a/test/ReachSets/unit_discretization.jl +++ b/test/ReachSets/unit_discretization.jl @@ -4,69 +4,76 @@ let # preventing global scope # =================================================================== A = sparse([1, 1, 2, 3, 4], [1, 2, 2, 4, 3], [1., 2., 3., 4., 5.], 4, 4) X0 = BallInf(zeros(4), 0.1) - #U = ZeroSet(size(A, 1)) - cont_sys_homog = ContinuousSystem(A, X0) + cont_sys_homog = IVP(LCS(A), X0) δ = 0.01 - # no bloating, do not use Pade approximation - discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="nobloating", pade_expm=false) + # no bloating, use default options + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="nobloating") @test inputdim(discr_sys_homog) == 0 # there is no input set! # no bloating, use Pade approximation - discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="nobloating", pade_expm=true) + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="nobloating", exp_method="pade") @test inputdim(discr_sys_homog) == 0 - # bloating, do not use Pade approximation - discr_sys_homog = discretize(cont_sys_homog, δ, pade_expm=false) + # bloating, default options + discr_sys_homog = discretize(cont_sys_homog, δ) @test inputdim(discr_sys_homog) == 0 # bloating, first order - discr_sys_homog = discretize(cont_sys_homog, δ, approx_model="firstorder") + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="firstorder") + @test inputdim(discr_sys_homog) == 0 + + # bloating, forward interpolation + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="forward") + @test inputdim(discr_sys_homog) == 0 + + # bloating, backward interpolation + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="backward") @test inputdim(discr_sys_homog) == 0 # =============================================================== # Discretization of a continuous-time system with constant input # =============================================================== - U = Ball2(ones(4), 0.5) - cont_sys = ContinuousSystem(A, X0, U) + U = ConstantInput(Ball2(ones(4), 0.5)) + B = Matrix{Float64}(1.0I, 4, 4) + cont_sys = IVP(CLCCS(A, B, nothing, U), X0) - # no bloating, do not use Pade approximation - discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=false) + # no bloating + discr_sys = discretize(cont_sys, δ, algorithm="nobloating") @test inputdim(discr_sys) == 4 - #@test length(discr_sys.U) == 1 # a constant input has no length method! - inputs = next_set(inputset(discr_sys)) # getter: inputset(discr_sys) = discr_sys.s.U + + inputs = next_set(inputset(discr_sys)) @test dim(inputs) == 4 @test isa(inputs, LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == ones(4) && inputs.X.radius == 0.5 # no bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=true) + discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade") @test inputdim(discr_sys) == 4 - # bloating, do not use Pade approximation - discr_sys = discretize(cont_sys, δ, pade_expm=false) + # bloating, use scaling and squaring method + discr_sys = discretize(cont_sys, δ, exp_method="base") @test inputdim(discr_sys) == 4 - #@test length(discr_sys.U) == 1 inputs = next_set(inputset(discr_sys)) @test dim(inputs) == 4 @test isa(inputs, MinkowskiSum) # bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, pade_expm=true) + discr_sys = discretize(cont_sys, δ, exp_method="pade") @test inputdim(discr_sys) == 4 - discr_sys = discretize(cont_sys, δ, approx_model="firstorder") + discr_sys = discretize(cont_sys, δ, algorithm="firstorder") @test inputdim(discr_sys) == 4 # =================================================================== # Discretization of a continuous-time system with time-varying input # =================================================================== Ui = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3] - cont_sys = ContinuousSystem(A, X0, Ui) + cont_sys = IVP(CLCCS(A, Matrix(1.0I, 4, 4), nothing, VaryingInput(Ui)), X0) - # no bloating, do not use Pade approximation - discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=false) + # no bloating + discr_sys = discretize(cont_sys, δ, algorithm="nobloating") Ui_d = inputset(discr_sys) @test length(Ui_d) == 3 for (i, inputs) in enumerate(discr_sys.s.U) @@ -86,10 +93,10 @@ let # preventing global scope end # no bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=true) + discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade") - # bloating, do not use Pade approximation - discr_sys = discretize(cont_sys, δ, pade_expm=false) + # bloating + discr_sys = discretize(cont_sys, δ, exp_method="base") Ui_d = inputset(discr_sys) @test length(Ui_d) == 3 for (i, inputs) in enumerate(discr_sys.s.U) @@ -112,9 +119,11 @@ let # preventing global scope end # bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, pade_expm=true) + discr_sys = discretize(cont_sys, δ, exp_method="pade") # lazy symmetric interval hull - discr_sys = discretize(cont_sys, δ, lazy_sih=true) - discr_sys = discretize(cont_sys, δ, lazy_sih=false) + discr_sys = discretize(cont_sys, δ, sih_method="lazy") + + # concrete symmetric interval hull + discr_sys = discretize(cont_sys, δ, sih_method="concrete") end diff --git a/test/Reachability/models/bouncing_ball.jl b/test/Reachability/models/bouncing_ball.jl index 659457a0..53a44392 100644 --- a/test/Reachability/models/bouncing_ball.jl +++ b/test/Reachability/models/bouncing_ball.jl @@ -16,7 +16,7 @@ function bouncing_ball() U = Singleton([1.0]) inv = HalfSpace([-1.0, 0.0], 0.0) # x >= 0 m1 = ConstrainedLinearControlContinuousSystem( - A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U) + A, Matrix(1.0I, size(B, 1), size(B, 1)), inv, ConstantInput(B*U)) # modes modes = [m1] diff --git a/test/Reachability/models/thermostat.jl b/test/Reachability/models/thermostat.jl index 9d1d34a8..9c64310d 100644 --- a/test/Reachability/models/thermostat.jl +++ b/test/Reachability/models/thermostat.jl @@ -13,21 +13,21 @@ function thermostat() # automaton structure automaton = LightAutomaton(2) + I1 = Matrix(1.0I, 1, 1) + # mode on A = hcat(-c_a) B = hcat(30.) U = Singleton([c_a]) inv = HalfSpace([1.0], 22.0) # x ≤ 22 - m_on = ConstrainedLinearControlContinuousSystem( - A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U) + m_on = ConstrainedLinearControlContinuousSystem(A, I1, inv, ConstantInput(B*U)) # mode off A = hcat(-c_a) B = hcat(0.0) U = Singleton([0.0]) inv = HalfSpace([-1.0], -18.0) # x ≥ 18 - m_off = ConstrainedLinearControlContinuousSystem( - A, Matrix{eltype(A)}(I, size(B, 1), size(B, 1)), inv, B*U) + m_off = ConstrainedLinearControlContinuousSystem(A, I1, inv, ConstantInput(B*U)) # modes modes = [m_on, m_off] diff --git a/test/Reachability/solve_continuous.jl b/test/Reachability/solve_continuous.jl index a429ccae..bddea11f 100644 --- a/test/Reachability/solve_continuous.jl +++ b/test/Reachability/solve_continuous.jl @@ -12,17 +12,17 @@ A = [ 0.0509836 0.168159 0.95246 0.33644 X0 = BallInf(ones(4), 0.1) # default options (computes all variables) -s = solve(InitialValueProblem(LinearContinuousSystem(A), X0), :T=>0.1) +s = solve(IVP(LCS(A), X0), :T=>0.1) @test dim(s.Xk[1].X) == 4 # two variables and custom partition -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1), op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4])) @test dim(s.Xk[1].X) == 4 # the default partition is used. -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1), op=BFFPSV18(:vars=>[1,3])) @test dim(s.Xk[1].X) == 2 @@ -30,7 +30,7 @@ s = solve(ContinuousSystem(A, X0), # =============================== # Test projection # =============================== -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1, :project_reachset=>false), op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4])) ps = project(s) @@ -45,14 +45,14 @@ A = [ 0.0509836 0.168159 0.95246 0.33644 0.38318 0.616014 0.518412 0.778765] # check that x1 + x2 <= 2 doesn't hold -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1, :mode=>"check", :property=>SafeStatesProperty(LinearConstraint([1., 1., 0., 0.], 2.))), op=BFFPSV18(:vars=>[1,2,3], :partition=>[1:2, 3:4])) @test s.violation == 1 # check that x1 - x2 <= 2 holds -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1, :mode=>"check", :property=>SafeStatesProperty(LinearConstraint([1., -1., 0., 0.], 2.))), op=BFFPSV18(:vars=>[1,2,3], :partition=>[1:2, 3:4])) @@ -63,42 +63,42 @@ s = solve(ContinuousSystem(A, X0), # =============================== # template directions (eg. :box, :oct, :octbox) -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1, :ε_proj=>1e-5, :set_type_proj=>HPolygon), op=BFFPSV18(:vars=>[1,3], :partition=>[1:4], :template_directions => :oct)) -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :lazy_sih=>true)) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :sih_method=>"lazy")) -s = solve(ContinuousSystem(A, X0), +s = solve(IVP(LCS(A), X0), Options(:T=>0.1), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :lazy_sih=>false)) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :sih_method=>"concrete")) -s = solve(ContinuousSystem(sparse(A), X0), +s = solve(IVP(LCS(sparse(A)), X0), Options(:T=>0.1), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :lazy_expm=>true, - :lazy_expm_discretize=>true, :assume_sparse=>false)) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :exp_method=>"lazy", + :assume_sparse=>false)) -s = solve(ContinuousSystem(sparse(A), X0), +s = solve(IVP(LCS(sparse(A)), X0), Options(:T=>0.1), - op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :lazy_expm=>true, - :lazy_expm_discretize=>true, :assume_sparse=>true)) + op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :exp_method=>"lazy", + :assume_sparse=>true)) # uses Interval for set_type_init and set_type_iter but Hyperrectangle for # set_type_proj -s = solve(ContinuousSystem(sparse(A), X0), +s = solve(IVP(LCS(sparse(A)), X0), Options(:T=>0.1), op=BFFPSV18(:vars=>[1,3], :partition=>[[i] for i in 1:4], - :lazy_expm=>true, :lazy_expm_discretize=>true, + :exp_method=>"lazy", :set_type=>Interval)) # =============================== # System with an odd dimension # =============================== A = randn(5, 5); X0 = BallInf(ones(5), 0.1) -s = solve(ContinuousSystem(A, X0), Options(:T=>0.1), +s = solve(IVP(LCS(sparse(A)), X0), Options(:T=>0.1), op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4, [5]], :block_types=> Dict{Type{<:LazySet}, AbstractVector{<:AbstractVector{Int}}}( HPolygon=>[1:2, 3:4], Interval=>[[5]]))) @@ -108,7 +108,7 @@ s = solve(ContinuousSystem(A, X0), Options(:T=>0.1), # =============================== A = [1 2 1.; 0 0. 1; -2 1 4] X0 = BallInf(ones(3), 0.1) -system = ContinuousSystem(A, X0) +system = IVP(LCS(sparse(A)), X0) s = solve(system, Options(:T=>0.1), op=BFFPSV18(:vars=>1:3, :partition=>[1:2, 3:3])) @@ -133,7 +133,7 @@ U = Interval(0.99, 1.01) × Interval(0.99, 1.01) X0 = BallInf(ones(4), 0.1) # dense identity matrix -E = Matrix{Float64}(I, size(A)) +E = Matrix(1.0I, size(A)) # instantiate system Δ = ConstrainedLinearControlContinuousSystem(A, E, nothing, ConstantInput(B*U)) diff --git a/test/Systems/unit_System.jl b/test/Systems/unit_System.jl index bf14e749..2a5a3367 100644 --- a/test/Systems/unit_System.jl +++ b/test/Systems/unit_System.jl @@ -4,7 +4,7 @@ let # preventing global scope # ====================================================== A = sparse([1, 1, 2, 3, 4], [1, 2, 2, 4, 3], [1., 2., 3., 4., 5.], 4, 4) X0 = BallInf(zeros(4), 0.1) - cont_sys_homog = ContinuousSystem(A, X0) + cont_sys_homog = IVP(LCS(A), X0) # Check if the input is constant @test inputdim(cont_sys_homog) == 0 @@ -17,8 +17,8 @@ let # preventing global scope # =================================================== # Testing continuous-time system with constant input # =================================================== - U = Ball2(ones(4), 0.5) - cont_sys = ContinuousSystem(A, X0, U) + U = ConstantInput(Ball2(ones(4), 0.5)) + cont_sys = IVP(CLCCS(A, Matrix(1.0I, 4, 4), nothing, U), X0) # check initial state @test cont_sys.x0.center ≈ zeros(4) && cont_sys.x0.radius ≈ 0.1 @@ -31,8 +31,8 @@ let # preventing global scope # ======================================================== # Testing continuous-time system with time-varying input # ======================================================== - Ui = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3] - cont_sys = ContinuousSystem(A, X0, Ui) + U = VaryingInput([Ball2(0.01*i*ones(4), i*0.2) for i in 1:3]) + cont_sys = cont_sys = IVP(CLCCS(A, Matrix(1.0I, 4, 4), nothing, U), X0) for (i, inputs) in enumerate(cont_sys.s.U) if i == 1 @test inputs.center ≈ 0.01*ones(4) @@ -50,7 +50,7 @@ let # preventing global scope # Testing discrete-time homogeneous system # ========================================= δ = 0.01 - discr_sys_homog = DiscreteSystem(A, X0) + discr_sys_homog = IVP(LDS(A), X0) # Check if the input is constant #@test isa(discr_sys_homog.s.U, ConstantInput) # no field U @@ -63,8 +63,8 @@ let # preventing global scope # ================================================= # Testing discrete-time system with constant input # ================================================= - U = Ball2(ones(4), 0.5) - discr_sys = DiscreteSystem(A, X0, U) + U = ConstantInput(Ball2(ones(4), 0.5)) + discr_sys = IVP(CLCDS(A, Matrix(1.0I, 4, 4), nothing, U), X0) # check initial state @test discr_sys.x0.center ≈ zeros(4) && discr_sys.x0.radius ≈ 0.1 @@ -77,8 +77,8 @@ let # preventing global scope # ===================================================== # Testing discrete-time system with time-varying input # ===================================================== - Ui = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3] - discr_sys = DiscreteSystem(A, X0, Ui) + U = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3] + discr_sys = IVP(CLCDS(A, Matrix(1.0I, 4, 4), nothing, U), X0) Uset = inputset(discr_sys.s) for (i, inputs) in enumerate(cont_sys.s.U) if i == 1 diff --git a/test/runtests.jl b/test/runtests.jl index adac2556..3086ec3d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,13 @@ using LazySets, Reachability, MathematicalSystems, LinearAlgebra, SparseArrays # fix namespace conflicts with MathematicalSystems -import LazySets.LinearMap +using LazySets: LinearMap + +# common aliases +const CLCCS = ConstrainedLinearControlContinuousSystem +const CLCDS = ConstrainedLinearControlDiscreteSystem +const LCS = LinearContinuousSystem +const LDS = LinearDiscreteSystem # compatibility between Julia versions include("../src/compat.jl")