From 3aec3a46b6c56693852d6ec46ca8425a324b7db1 Mon Sep 17 00:00:00 2001 From: mforets Date: Sat, 2 Mar 2019 00:20:02 -0300 Subject: [PATCH 01/24] revise discretize --- src/ReachSets/discretize.jl | 246 +++++++++++++++++++++++++----------- 1 file changed, 169 insertions(+), 77 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 23223202..b71e4b57 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -1,97 +1,185 @@ """ - discretize(cont_sys, δ; [approx_model], [pade_expm], [lazy_expm], [lazy_sih]) + discretize(𝑆, δ; [approx_model], [expm_method], [sih_method]) -Discretize a continuous system of ODEs with nondeterministic inputs. +Apply an approximation model to `S` obtaining a discrete initial value problem. -## Input +### Input -- `cont_sys` -- continuous system +- `𝑆` -- initial value problem for a continuous affine ODE with + nondeterministic inputs - `δ` -- step size - `approx_model` -- the method to compute the approximation model for the - discretization, among: + 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 - (use for discrete-time reachability) -- `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) +- `expm_method` -- (optional, default: `base`) the method used to take the matrix + exponential of the coefficients 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, as implemented in `Expokit` + - `lazy` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` -A discrete system. +- `sih_method` -- (optional, default: `lazy`) 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 of a discrete system. + +### Algorithm + +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 nondeterministic inputs +and ``\\mathcal{X}_0`` is the set of initial states. Recall that the system +`𝑆` is called homogeneous whenever `U` is the empty set. -## Notes +Given a step size ``δ``, this function computes a set, `Ω0`, 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, δ]``. -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"`. +The initial value problem returned by this function consists of the set `Ω0` +together with the coefficients matrix ``ϕ = e^{Aδ}`` and a transformed +set of inputs if `U` is non-empty. + +In the literature, the method to obtain `Ω0` is called the *approximation model* +and different alternatives have been proposed. See the argument `approx_model` +for available options. For the reference to the original papers, see the docstring +of each method. + +The transformation described above 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. In particular, if you don't want to consider +bloating, i.e. discrete-time reachability, use `approx_model="nobloating"`. + +Several methods to compute the matrix exponential are availabe. Use `expm_method` +to select one. For very large systems (~10000×10000), computing the full matrix +exponential is very expensive hence it is preferable to compute the action +of the matrix exponential over vectors when needed. Use the option +`expm_method="lazy"` for this. """ -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} + expm_method::String="base", + sih_method::String="lazy") + + if !isaffine(S) + throw(ArgumentError("`discretize` is only implemented for affine ODEs")) + end if approx_model in ["forward", "backward"] - return discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm, - lazy_expm, lazy_sih) + return discr_bloat_interpolation(𝑆, δ, approx_model, expm_method, sih_method) elseif approx_model == "firstorder" - return discr_bloat_firstorder(cont_sys, δ) + return _discretize_first_order(𝑆, δ, expm_method) elseif approx_model == "nobloating" - return discr_no_bloat(cont_sys, δ, pade_expm, lazy_expm) + return discr_no_bloat(𝑆, δ, expm_method) else - error("The approximation model is invalid") + throw(ArgumentError("the approximation model $approx_model is unknown")) + end +end + +function exp_Aδ(A::AbstractMatrix, δ::Float64, expm_method="base") + if expm_method == "base" + return expmat(Matrix(A*δ)) + elseif expm_method == "lazy" + return SparseMatrixExp(A*δ) + elseif expm_method == "pade" + return padm(A*δ) + else + throw(ArgumentError("the exponentiation method $expm_method is unknown")) end end """ - bloat_firstorder(cont_sys, δ) + _discretize_first_order(𝑆, δ, [p], [expm_method]) -Compute bloating factors using first order approximation. +Apply a first order approximation model to `S` obtaining a discrete initial value problem. -## Input +### Input -- `cont_sys` -- a continuous affine system -- `δ` -- step size +- `𝑆` -- initial value problem for a continuous affine ODE with nondeterministic inputs +- `δ` -- step size +- `p` -- (optional, default: `Inf`) parameter in the considered norm -## Notes +### Output -In this algorithm, the infinity norm is used. -See also: `discr_bloat_interpolation` for more accurate (less conservative) -bounds. +The initial value problem for a discrete system. -## Algorithm +### Algorithm + +Let us define some notation. 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 nondeterministic inputs and ``\\mathcal{X}_0`` is the set +of initial states. + +Let ``R_{X0} = \\max_{x∈ X0} ||x||`` and `D_{X0} = \\max_{x, y ∈ X0} ||x-y||`` +and ``R_{V} = \\max_{u ∈ U} ||u||``. + +It is proved [Lemma 1, 1] that the set `Ω0` defined as + +```math +Ω0 = CH(\\mathcal{X}_0, e^{δA}\\mathcal{X}_0 ⊕ δU ⊕ αB_p) +``` +where ``α = (e^{δ||A||} - 1 - δ||A||)*R_{X0} + R_{U} / ||A||)`` and ``B_p`` denotes +the unit ball for the considered norm, is such that the set of states reachable +by ``S`` in the time interval ``[0, δ]``, that we denote ``R_{[0,δ]}(X0)``, +is included in `Ω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(Ω0, R_{[0,δ]}(X0)) ≤ \\frac{1}{4}(e^{δ||A||} - 1) D_{X0} + 2α. +``` -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 +### Notes + +In this implementation, the infinity norm is used by default. To use other norms +substitute `BallInf` with the ball in the appropriate norm. However, note that +not all norms are supported; see the documentation of `?norm` in `LazySets` for +details. + +[1] Le Guernic, C., & Girard, A., 2010, *Reachability analysis of linear systems using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.* + +### Notes + +See also [`discr_bloat_interpolation`](@ref) for an alternative algorithm that +uses less conservative bounds. """ -function discr_bloat_firstorder(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, - δ::Float64) +function _discretize_first_order(𝑆::InitialValueProblem, δ::Float64, + p::Float64=Inf, + expm_method::String="base") + + # unwrap coeffs matrix and initial states + A, X0 = 𝑆.s.A, 𝑃.x0 + + # system size; A is assumed square + n = size(A, 1) - 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 + # compute exp(A*δ) + ϕ = exp_Aδ(A, δ, expm_method) + + if islinear(𝑆) # inputdim(𝑆) == 0 α = (exp(δ*Anorm) - 1. - δ*Anorm) * RX0 - Ω0 = CH(X0, ϕ * X0 + Ball2(zeros(size(ϕ, 1)), α)) - return DiscreteSystem(ϕ, Ω0) + □ = Ballp(p, zeros(n), α) + Ω0 = CH(X0, ϕ * X0 + □) + return InitialValueProblem(LinearDiscreteSystem(ϕ), Ω0)) # @system x' = ϕ*x, x(0) ∈ Ω0 + # ---- TODO ---- seguir aca ---- else # affine case; TODO: unify Constant and Varying input branches? Uset = inputset(cont_sys) @@ -115,8 +203,6 @@ function discr_bloat_firstorder(cont_sys::InitialValueProblem{<:AbstractContinuo return DiscreteSystem(ϕ, Ω0, discr_U) end end - - end """ @@ -157,26 +243,23 @@ function discr_no_bloat(cont_sys::InitialValueProblem{<:AbstractContinuousSystem pade_expm::Bool, lazy_expm::Bool) + # unrwap coefficients matrix and initial states 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 + + # compute matrix ϕ = exp(Aδ) + ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) # early return for homogeneous systems if cont_sys isa IVP{<:LinearContinuousSystem} Ω0 = X0 return DiscreteSystem(ϕ, Ω0) end + U = inputset(cont_sys) inputs = next_set(U, 1) + n = size(A, 1) + # compute matrix to transform the inputs if lazy_expm P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); @@ -237,26 +320,19 @@ The matrix `P` is such that: `ϕAabs = P[1:n, 1:n]`, """ function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64, - approx_model::String, - pade_expm::Bool, - lazy_expm::Bool, - lazy_sih::Bool) + approx_model::String="forward", + pade_expm::Bool=false, + lazy_expm::Bool=false, + lazy_sih::Bool=true) sih = lazy_sih ? SymmetricIntervalHull : symmetric_interval_hull + # unrwap coefficients matrix and initial states A, X0 = cont_sys.s.A, cont_sys.x0 n = size(A, 1) # compute matrix ϕ = exp(Aδ) - if lazy_expm - ϕ = SparseMatrixExp(A*δ) - else - if pade_expm - ϕ = padm(A*δ) - else - ϕ = expmat(Matrix(A*δ)) - end - end + ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) # early return for homogeneous systems if cont_sys isa IVP{<:LinearContinuousSystem} @@ -308,3 +384,19 @@ function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:AbstractConti return DiscreteSystem(ϕ, Ω0, discretized_U) end end + +function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:LinearContinuousSystem}, + δ::Float64, + approx_model::String="forward", + pade_expm::Bool=false, + lazy_expm::Bool=false, + lazy_sih::Bool=true) + + # unwrap coefficients matrix and initial states + A, X0 = cont_sys.s.A, cont_sys.x0 + + # compute matrix ϕ = exp(Aδ) + ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) + + +end From bb9deff7e2a963da72e401497aa9888e7a586f75 Mon Sep 17 00:00:00 2001 From: mforets Date: Sun, 3 Mar 2019 23:43:17 -0300 Subject: [PATCH 02/24] update branch (wip) --- src/ReachSets/discretize.jl | 540 ++++++++++++++++++++++-------------- 1 file changed, 337 insertions(+), 203 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index b71e4b57..c5503de0 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -1,39 +1,45 @@ +const LDS = LinearDiscreteSystem +const CLCDS = ConstrainedLinearControlDiscreteSystem + +@inline I(T, n) = Matrix{eltype(A)}(I, n, n) + """ - discretize(𝑆, δ; [approx_model], [expm_method], [sih_method]) + discretize(𝑆, δ; [approximation], [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 - nondeterministic inputs -- `δ` -- step size -- `approx_model` -- the method 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 - -- `expm_method` -- (optional, default: `base`) the method used to take the matrix - exponential of the coefficients 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, as implemented in `Expokit` - - `lazy` -- compute a wrapper type around the matrix exponential, i.e. using - the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` - -- `sih_method` -- (optional, default: `lazy`) the method used to take the +- `𝑆` -- initial value problem for a continuous affine ODE with + non-deterministic inputs +- `δ` -- step size +- `approximation` -- the method 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: `"lazy"`) 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 + - `"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` ### Output @@ -42,77 +48,218 @@ The initial value problem of a discrete system. ### Algorithm 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 nondeterministic inputs +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. -Given a step size ``δ``, this function computes a set, `Ω0`, that guarantees to +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, δ]``. -The initial value problem returned by this function consists of the set `Ω0` -together with the coefficients matrix ``ϕ = e^{Aδ}`` and a transformed +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. -In the literature, the method to obtain `Ω0` is called the *approximation model* -and different alternatives have been proposed. See the argument `approx_model` +In the literature, the method to obtain `Ω₀` is called the *approximation model* +and different alternatives have been proposed. See the argument `approximation` for available options. For the reference to the original papers, see the docstring of each method. -The transformation described above 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. In particular, if you don't want to consider -bloating, i.e. discrete-time reachability, use `approx_model="nobloating"`. +In the dense-time case, the transformation described is such that the trajectories +of the given continuous system are included in the computed flowpipe of the +discretized system. + +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 +`approximation="nobloating"` for this setting. -Several methods to compute the matrix exponential are availabe. Use `expm_method` +Several methods to compute the matrix exponential are availabe. Use `exp_method` to select one. For very large systems (~10000×10000), computing the full matrix exponential is very expensive hence it is preferable to compute the action of the matrix exponential over vectors when needed. Use the option -`expm_method="lazy"` for this. +`exp_method="lazy"` for this. """ function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; - approx_model::String="forward", - expm_method::String="base", + approximation::String="forward", + exp_method::String="base", sih_method::String="lazy") - if !isaffine(S) - throw(ArgumentError("`discretize` is only implemented for affine ODEs")) - end - - if approx_model in ["forward", "backward"] - return discr_bloat_interpolation(𝑆, δ, approx_model, expm_method, sih_method) - elseif approx_model == "firstorder" - return _discretize_first_order(𝑆, δ, expm_method) - elseif approx_model == "nobloating" - return discr_no_bloat(𝑆, δ, expm_method) + if approximation in ["forward", "backward"] + return _discretize_interpolation(𝑆, δ, approximation=approximation, + exp_method=exp_method, sih_method=sih_method) + elseif approximation == "firstorder" + return _discretize_firstorder(𝑆, δ, exp_method=exp_method) + elseif approximation == "nobloating" + return _discretize_nobloating(𝑆, δ, exp_method=exp_method) else - throw(ArgumentError("the approximation model $approx_model is unknown")) + throw(ArgumentError("the approximation model $approximation is unknown")) end end -function exp_Aδ(A::AbstractMatrix, δ::Float64, expm_method="base") - if expm_method == "base" +""" + exp_Aδ(A::AbstractMatrix, δ::Float64; [exp_method]) + +Compute the matrix exponential ``e^{Aδ}``. + +### 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, as implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` + +### Output + +A matrix. +""" +function exp_Aδ(A::AbstractMatrix{Float64}, δ::Float64; exp_method="base") + if exp_method == "base" return expmat(Matrix(A*δ)) - elseif expm_method == "lazy" + elseif exp_method == "lazy" return SparseMatrixExp(A*δ) - elseif expm_method == "pade" + elseif exp_method == "pade" return padm(A*δ) else - throw(ArgumentError("the exponentiation method $expm_method is unknown")) + throw(ArgumentError("the exponentiation method $exp_method is unknown")) + end +end + +""" + ϕ₁(A, δ; [exp_method]) + +TODO: Add doctring + +### 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, as implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` + +### Output + +A matrix. +""" +function ϕ₁(A, δ; exp_method="base") + n = size(A, 1) + if exp_method == "base" + P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*n)])) + ϕ₁_Aδ = P[1:n, (n+1):2*n] + + elseif exp_method == "lazy" + P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*n)]) + ϕ₁_Aδ = sparse(get_columns(P, (n+1):2*n)[1:n, :]) + + elseif exp_method == "pade" + P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*n)]) + ϕ₁_Aδ = P[1:n, (n+1):2*n] + + else + throw(ArgumentError("the exponentiation method $exp_method is unknown")) + end + + return ϕ₁_Aδ +end + +""" + ϕ₂(A, δ; [exp_method]) + +TODO: Add doctring + +### 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, as implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` + +### Output + +A matrix. +""" +function ϕ₂(A, δ; exp_method="base") + n = size(A, 1) + if exp_method == "base" + P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*n)])) + ϕ₂_Aδ = P[1:n, (2*n+1):3*n] + + elseif exp_method == "lazy" + P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*n)]) + ϕ₂_Aδ = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) + + elseif exp_method == "pade" + P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); + spzeros(n, 2*n) sparse(δ*I, n, n); + spzeros(n, 3*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_first_order(𝑆, δ, [p], [expm_method]) + _discretize_firstorder(𝑆, δ; [p], [exp_method]) -Apply a first order approximation model to `S` obtaining a discrete initial value problem. +Apply a first-order approximation model to `S` obtaining a discrete initial value problem. ### Input -- `𝑆` -- initial value problem for a continuous affine ODE with nondeterministic inputs -- `δ` -- step size -- `p` -- (optional, default: `Inf`) parameter in the considered norm +- `𝑆` -- 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, as implemented in `Expokit` + - `lazy` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` ### Output @@ -122,25 +269,33 @@ The initial value problem for a discrete system. Let us define some notation. 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 nondeterministic inputs and ``\\mathcal{X}_0`` is the set +where `U` is the set of non-deterministic inputs and ``\\mathcal{X}_0`` is the set of initial states. -Let ``R_{X0} = \\max_{x∈ X0} ||x||`` and `D_{X0} = \\max_{x, y ∈ X0} ||x-y||`` -and ``R_{V} = \\max_{u ∈ U} ||u||``. +Let ``R_{\\mathcal{X}_0} = \\max_{x ∈ \\mathcal{X}_0} ‖x‖``, +`D_{\\mathcal{X}_0} = \\max_{x, y ∈ \\mathcal{X}_0} ‖x-y‖`` and +``R_{V} = \\max_{u ∈ U} ‖u‖``. + +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 norm. -It is proved [Lemma 1, 1] that the set `Ω0` defined as +It is proved in [Lemma 1, 1] that the set of states reachable by ``S`` in the time +interval ``[0, δ]``, that we denote ``R_{[0,δ]}(\\mathcal{X}_0)``, +is included in ``Ω₀``: ```math -Ω0 = CH(\\mathcal{X}_0, e^{δA}\\mathcal{X}_0 ⊕ δU ⊕ αB_p) +R_{[0,δ]}(\\mathcal{X}_0) ⊆ Ω₀. ``` -where ``α = (e^{δ||A||} - 1 - δ||A||)*R_{X0} + R_{U} / ||A||)`` and ``B_p`` denotes -the unit ball for the considered norm, is such that the set of states reachable -by ``S`` in the time interval ``[0, δ]``, that we denote ``R_{[0,δ]}(X0)``, -is included in `Ω0`. Moreover, if `d_H(A, B)` denotes the Hausdorff distance -between the sets ``A`` and ``B`` in ``\\mathbb{R}^n``, then + +Moreover, if `d_H(A, B)` denotes the Hausdorff distance between the sets ``A`` +and ``B`` in ``\\mathbb{R}^n``, then ```math -d_H(Ω0, R_{[0,δ]}(X0)) ≤ \\frac{1}{4}(e^{δ||A||} - 1) D_{X0} + 2α. +d_H(Ω₀, R_{[0,δ]}(\\mathcal{X}_0)) ≤ \\frac{1}{4}(e^{δ ‖A‖} - 1) D_{\\mathcal{X}_0} + 2α. ``` ### Notes @@ -150,75 +305,87 @@ substitute `BallInf` with the ball in the appropriate norm. However, note that not all norms are supported; see the documentation of `?norm` in `LazySets` for details. -[1] Le Guernic, C., & Girard, A., 2010, *Reachability analysis of linear systems -using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.* - -### Notes - See also [`discr_bloat_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_first_order(𝑆::InitialValueProblem, δ::Float64, - p::Float64=Inf, - expm_method::String="base") +function _discretize_firstorder(𝑆::InitialValueProblem, + δ::Float64; + p::Float64=Inf, + exp_method::String="base") - # unwrap coeffs matrix and initial states - A, X0 = 𝑆.s.A, 𝑃.x0 + # 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), Inf) - RX0 = norm(X0, Inf) + Anorm = norm(Matrix(A), p) + RX0 = norm(X0, p) # compute exp(A*δ) - ϕ = exp_Aδ(A, δ, expm_method) + ϕ = exp_Aδ(A, δ, exp_method) if islinear(𝑆) # inputdim(𝑆) == 0 α = (exp(δ*Anorm) - 1. - δ*Anorm) * RX0 □ = Ballp(p, zeros(n), α) - Ω0 = CH(X0, ϕ * X0 + □) - return InitialValueProblem(LinearDiscreteSystem(ϕ), Ω0)) # @system x' = ϕ*x, x(0) ∈ Ω0 - # ---- TODO ---- seguir aca ---- - else - # affine case; TODO: unify Constant and Varying input branches? - Uset = inputset(cont_sys) + Ω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) + α = (exp(δ*Anorm) - 1.0 - δ*Anorm)*(RX0 + RU/Anorm) + β = (exp(δ*Anorm) - 1.0 - δ*Anorm)*RU/Anorm + □α = Ballp(p, zeros(n), α) + □β = Ballp(p, zeros(n), β) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ * U + □α) + Ud = map(u -> δ*u ⊕ □β, U) + return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) + elseif Uset isa VaryingInput - discr_U = Vector{LazySet}(undef, length(Uset)) + Ud = Vector{LazySet}(undef, length(Uset)) # TODO: concrete type of 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) + α = (exp(δ*Anorm) - 1.0 - δ*Anorm)*(RX0 + RU/Anorm) + β = (exp(δ*Anorm) - 1.0 - δ*Anorm)*RU/Anorm + □α = Ballp(p, zeros(n), α) + □β = Ballp(p, zeros(n), β) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ * Ui ⊕ □α) + Ud[i] = δ * Ui ⊕ □β end - return DiscreteSystem(ϕ, Ω0, discr_U) + Ud = VaryingInput(Ud) + return IVP(CLCDS(ϕ, I(typeof(ϕ), n), nothing, Ud), Ω0) 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 -- `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: + + - `"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, as implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` ## Output @@ -238,46 +405,27 @@ Verification of Hybrid Systems.* 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. """ -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") - # unrwap coefficients matrix and initial states - A, X0 = cont_sys.s.A, cont_sys.x0 + # unrwap coefficient matrix and initial states + A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) # early return for homogeneous systems - if cont_sys isa IVP{<:LinearContinuousSystem} + if islinear(𝑆) Ω0 = X0 - return DiscreteSystem(ϕ, Ω0) + return IVP(LDS(ϕ), Ω0) end - U = inputset(cont_sys) + U = inputset(𝑆) inputs = next_set(U, 1) - n = size(A, 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 + Phi1Adelta = ϕ₁(A, δ, exp_method) discretized_U = Phi1Adelta * inputs @@ -292,20 +440,34 @@ function discr_no_bloat(cont_sys::InitialValueProblem{<:AbstractContinuousSystem end """ - discr_bloat_interpolation(cont_sys, δ, approx_model, pade_expm, lazy_expm) + _discretize_interpolation(𝑆, δ, [approximation], [exp_method], [sih_method]) Compute bloating factors using forward or backward interpolation. ## Input -- `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) +- `cs` -- a continuous system +- `δ` -- step size +- `approximation` -- choose 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: + + - `"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, as implemented in `Expokit` + - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using + the lazy implementation `SparseMatrixExp` from `LazySets` and + evaluation of the action of the matrix exponential using the + `expmv` implementation in `Expokit` + +- `sih_method` -- (optional, default: `"lazy"`) 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 ## Algorithm @@ -318,62 +480,50 @@ be obtained directly, as a function of the inverse of A and `e^{At} - I`. 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="forward", - pade_expm::Bool=false, - lazy_expm::Bool=false, - lazy_sih::Bool=true) - - sih = lazy_sih ? SymmetricIntervalHull : symmetric_interval_hull +function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, + δ::Float64; + approximation::String="forward", + exp_method::String="base", + sih_method::String="lazy") + + if sih_method == "lazy" + sih = SymmetricIntervalHull + elseif sih_method == "concrete" + sih = symmetric_interval_hull + else + throw(ArgumentError("the method $sih_method is unknown")) + end - # unrwap coefficients matrix and initial states - A, X0 = cont_sys.s.A, cont_sys.x0 - n = size(A, 1) + # unrwap coefficient matrix and initial states + A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) # early return for homogeneous systems - if cont_sys isa IVP{<:LinearContinuousSystem} - Ω0 = CH(X0, ϕ * X0) - return DiscreteSystem(ϕ, Ω0) + if islinear(𝑆) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ E) + return IVP(LDS(ϕ), Ω0) end - U = inputset(cont_sys) + U = inputset(𝑆) inputs = next_set(U, 1) # 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, :]) - 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] - end + Phi2Aabs = ϕ₂_Aδ(abs.(A), δ, exp_method=exp_method) if isa(inputs, ZeroSet) - if approx_model == "forward" || approx_model == "backward" - Ω0 = CH(X0, ϕ * X0 + δ * inputs) + if approximation == "forward" || approximation == "backward" + Ω0 = ConvexHull(X0, ϕ * X0 + δ * inputs) end else EPsi = sih(Phi2Aabs * sih(A * inputs)) discretized_U = δ * inputs + EPsi - if approx_model == "forward" + if approximation == "forward" EOmegaPlus = sih(Phi2Aabs * sih((A * A) * X0)) - Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaPlus) - elseif approx_model == "backward" + Ω0 = ConvexHull(X0, ϕ * X0 + discretized_U + EOmegaPlus) + elseif approximation == "backward" EOmegaMinus = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) - Ω0 = CH(X0, ϕ * X0 + discretized_U + EOmegaMinus) + Ω0 = ConvexHull(X0, ϕ * X0 + discretized_U + EOmegaMinus) end end @@ -384,19 +534,3 @@ function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:AbstractConti return DiscreteSystem(ϕ, Ω0, discretized_U) end end - -function discr_bloat_interpolation(cont_sys::InitialValueProblem{<:LinearContinuousSystem}, - δ::Float64, - approx_model::String="forward", - pade_expm::Bool=false, - lazy_expm::Bool=false, - lazy_sih::Bool=true) - - # unwrap coefficients matrix and initial states - A, X0 = cont_sys.s.A, cont_sys.x0 - - # compute matrix ϕ = exp(Aδ) - ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) - - -end From dbeefef0a07d1e2401eb973e91d49a82ee8a2d96 Mon Sep 17 00:00:00 2001 From: mforets Date: Mon, 4 Mar 2019 16:28:07 -0300 Subject: [PATCH 03/24] updates in discretizes --- src/ReachSets/discretize.jl | 413 +++++++++++++++++++++++------------- 1 file changed, 269 insertions(+), 144 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index c5503de0..58a5f987 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -65,7 +65,7 @@ and different alternatives have been proposed. See the argument `approximation` for available options. For the reference to the original papers, see the docstring of each method. -In the dense-time case, the transformation described is such that the trajectories +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. @@ -74,10 +74,10 @@ input is assumed to remain constant between sampled times. Use the option `approximation="nobloating"` for this setting. Several methods to compute the matrix exponential are availabe. Use `exp_method` -to select one. For very large systems (~10000×10000), computing the full matrix -exponential is very expensive hence it is preferable to compute the action -of the matrix exponential over vectors when needed. Use the option -`exp_method="lazy"` for this. +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(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; @@ -112,11 +112,11 @@ Compute the matrix exponential ``e^{Aδ}``. - `"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, as implemented in `Expokit` + sparse matrices, implemented in `Expokit` - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` ### Output @@ -135,9 +135,14 @@ function exp_Aδ(A::AbstractMatrix{Float64}, δ::Float64; exp_method="base") end """ - ϕ₁(A, δ; [exp_method]) + Φ₁(A, δ; [exp_method]) + +Compute the series -TODO: Add doctring +```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 @@ -149,47 +154,75 @@ TODO: Add doctring - `"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, as implemented in `Expokit` + sparse matrices, implemented in `Expokit` - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` - + 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 ϕ₁(A, δ; exp_method="base") +function Φ₁(A, δ; exp_method="base") n = size(A, 1) if exp_method == "base" P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)])) - ϕ₁_Aδ = P[1:n, (n+1):2*n] + Φ₁_Aδ = P[1:n, (n+1):2*n] elseif exp_method == "lazy" P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)]) - ϕ₁_Aδ = sparse(get_columns(P, (n+1):2*n)[1:n, :]) + Φ₁_Aδ = sparse(get_columns(P, (n+1):2*n)[1:n, :]) elseif exp_method == "pade" P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)]) - ϕ₁_Aδ = P[1:n, (n+1):2*n] + Φ₁_Aδ = P[1:n, (n+1):2*n] else - throw(ArgumentError("the exponentiation method $exp_method is unknown")) + throw(ArgumentError("the exponentiation method $exp_method is unknown")) end - return ϕ₁_Aδ + return Φ₁_Aδ end """ - ϕ₂(A, δ; [exp_method]) + Φ₂(A, δ; [exp_method]) -TODO: Add doctring +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 @@ -201,47 +234,72 @@ TODO: Add doctring - `"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, as implemented in `Expokit` + sparse matrices, implemented in `Expokit` - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` + 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") +function Φ₂(A, δ; exp_method="base") n = size(A, 1) if exp_method == "base" P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)])) - ϕ₂_Aδ = P[1:n, (2*n+1):3*n] + Φ₂_Aδ = P[1:n, (2*n+1):3*n] elseif exp_method == "lazy" P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)]) - ϕ₂_Aδ = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) + Φ₂_Aδ = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) elseif exp_method == "pade" P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); spzeros(n, 2*n) sparse(δ*I, n, n); spzeros(n, 3*n)]) - ϕ₂_Aδ = P[1:n, (2*n+1):3*n] + Φ₂_Aδ = P[1:n, (2*n+1):3*n] else throw(ArgumentError("the exponentiation method $exp_method is unknown")) end - return ϕ₂_Aδ + return Φ₂_Aδ end """ _discretize_firstorder(𝑆, δ; [p], [exp_method]) -Apply a first-order approximation model to `S` obtaining a discrete initial value problem. +Apply a first-order approximation model to `S` obtaining a discrete initial +value problem. ### Input @@ -249,63 +307,75 @@ Apply a first-order approximation model to `S` obtaining a discrete initial valu 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, as implemented in `Expokit` - - `lazy` -- compute a wrapper type around the matrix exponential, i.e. using - the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` +- `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. +The initial value problem for a discrete system. In particular: + +- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` + is returned, +- otherwise a constrained linear discrete systen is returned, + `ConstrainedLinearControlDiscreteSystem`. ### Algorithm -Let us define some notation. 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. +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. -Let ``R_{\\mathcal{X}_0} = \\max_{x ∈ \\mathcal{X}_0} ‖x‖``, +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_{V} = \\max_{u ∈ U} ‖u‖``. +``R_{V} = \\max_{u ∈ U} ‖u‖``. If only the support functions of ``\\mathcal{X}_0`` +and ``V`` 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 norm. +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, δ]``, that we denote ``R_{[0,δ]}(\\mathcal{X}_0)``, +interval ``[0, δ]``, that 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`` +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. ### Notes -In this implementation, the infinity norm is used by default. To use other norms -substitute `BallInf` with the ball in the appropriate norm. However, note that -not all norms are supported; see the documentation of `?norm` in `LazySets` for -details. +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 [`discr_bloat_interpolation`](@ref) for an alternative algorithm that +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 @@ -323,39 +393,55 @@ function _discretize_firstorder(𝑆::InitialValueProblem, 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_Aδ(A, δ, exp_method=exp_method) if islinear(𝑆) # inputdim(𝑆) == 0 - α = (exp(δ*Anorm) - 1. - δ*Anorm) * RX0 - □ = Ballp(p, zeros(n), α) - Ω0 = ConvexHull(X0, ϕ * X0 ⊕ □) + α = κ * 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.0 - δ*Anorm)*(RX0 + RU/Anorm) - β = (exp(δ*Anorm) - 1.0 - δ*Anorm)*RU/Anorm + 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), β) - Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ * U + □α) Ud = map(u -> δ*u ⊕ □β, U) return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) elseif Uset isa VaryingInput - Ud = Vector{LazySet}(undef, length(Uset)) # TODO: concrete type of Uset + Ud = Vector{LazySet}(undef, length(Uset)) # TODO: use concrete type of Uset? how? for (i, Ui) in enumerate(Uset) RU = norm(Ui, p) - α = (exp(δ*Anorm) - 1.0 - δ*Anorm)*(RX0 + RU/Anorm) - β = (exp(δ*Anorm) - 1.0 - δ*Anorm)*RU/Anorm - □α = Ballp(p, zeros(n), α) + + # 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), β) - Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ * Ui ⊕ □α) - Ud[i] = δ * Ui ⊕ □β + Ud[i] = δ*Ui ⊕ □β end Ud = VaryingInput(Ud) return IVP(CLCDS(ϕ, I(typeof(ϕ), n), nothing, Ud), Ω0) @@ -371,39 +457,54 @@ end Discretize a continuous system without bloating of the initial states, suitable for discrete-time reachability. -## Input +### Input - `𝑆` -- a continuous system - `δ` -- step size - `exp_method` -- (optional, default: `"base"`) the method used to take the matrix - exponential of the coefficient matrix, choose among: + 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, as implemented in `Expokit` + sparse matrices, implemented in `Expokit` - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` + the evaluation of the action of the matrix exponential using the + `expmv` implementation from `Expokit` -## Output +### Output -A discrete system. +The initial value problem for a discrete system. In particular: -## Algorithm +- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` + is returned, +- otherwise a constrained linear discrete systen is returned, + `ConstrainedLinearControlDiscreteSystem`. + +### Algorithm -The transformation implemented here is the following: +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. -- `A -> Phi := exp(A*delta)` -- `U -> V := M*U` -- `X0 -> X0hat := X0` +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. -where `M` corresponds to `Phi1(A, delta)` in Eq. (8) of *SpaceEx: Scalable -Verification of Hybrid Systems.* +The transformations are: -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. +- ``Φ ← \\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 _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; @@ -413,54 +514,46 @@ function _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousS A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) - ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) + ϕ = exp_Aδ(A, δ, exp_method=exp_method) + + # initial states remain unchanged + Ω0 = copy(X0) # early return for homogeneous systems if islinear(𝑆) - Ω0 = X0 return IVP(LDS(ϕ), Ω0) end - U = inputset(𝑆) - inputs = next_set(U, 1) - # compute matrix to transform the inputs - Phi1Adelta = ϕ₁(A, δ, exp_method) - - discretized_U = Phi1Adelta * inputs - - Ω0 = X0 + Phi1Adelta = Φ₁_Aδ(A, δ, exp_method=exp_method) + U = inputset(𝑆) + Ud = map(ui -> Phi1Adelta * ui, U) - 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(ϕ, I(typeof(A), n), nothing, Ud), Ω0) end """ - _discretize_interpolation(𝑆, δ, [approximation], [exp_method], [sih_method]) + _discretize_interpolation(𝑆, δ; [approximation], [exp_method], [sih_method]) Compute bloating factors using forward or backward interpolation. -## Input +### Input -- `cs` -- a continuous system +- `𝑆` -- a continuous system - `δ` -- step size - `approximation` -- choose 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: - - `"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, as implemented in `Expokit` - - `"lazy"` -- compute a wrapper type around the matrix exponential, i.e. using - the lazy implementation `SparseMatrixExp` from `LazySets` and - evaluation of the action of the matrix exponential using the - `expmv` implementation in `Expokit` + - `"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: `"lazy"`) the method used to take the symmetric interval hull operation, choose among: @@ -469,16 +562,44 @@ Compute bloating factors using forward or backward interpolation. - `"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 linear discrete system, `LinearDiscreteSystem` + is returned, +- otherwise a constrained linear discrete systen is returned, + `ConstrainedLinearControlDiscreteSystem`. + ## Algorithm -See Frehse et al., CAV'11, *SpaceEx: Scalable Verification of Hybrid Systems*, -Lemma 3. +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. + +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. + +The transformations are: -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`. +- ``Φ ← \\exp^{Aδ}``, +- ``Ω₀ ← ConvexHull(\\mathcal{X}_0, Φ\\mathcal{X}_0 ⊕ δU(0) ⊕ Eψ(U(0), δ) ⊕ E^+(\\mathcal{X}_0, δ))``, +- ``V ← δU(k) ⊕ Eψ(U(k), δ)``. -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]`. +Here we allow ``U`` to be a sequence of time varying non-deterministic input sets. + +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; @@ -498,39 +619,43 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) - ϕ = exp_Aδ(A, δ, lazy_expm, pade_expm) + ϕ = exp_Aδ(A, δ, exp_method=exp_method) + + # compute the transformation matrix to bloat the initial states + Phi2Aabs = Φ₂_Aδ(abs.(A), δ, exp_method=exp_method) + + if approximation == "forward" + Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use Eplus + elseif approximation == "backward" + Einit = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) # use Eminus + else + throw(ArgumentError("the method $approximation is unknown")) + end # early return for homogeneous systems if islinear(𝑆) - Ω0 = ConvexHull(X0, ϕ * X0 ⊕ E) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ Einit) return IVP(LDS(ϕ), Ω0) end + U = inputset(𝑆) - inputs = next_set(U, 1) + U0 = next_set(U, 1) - # compute the transformation matrix to bloat the initial states - Phi2Aabs = ϕ₂_Aδ(abs.(A), δ, exp_method=exp_method) + Eψ0 = sih(Phi2Aabs * sih(A * U0)) + Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ*U0 ⊕ Eψ0 ⊕ Einit) - if isa(inputs, ZeroSet) - if approximation == "forward" || approximation == "backward" - Ω0 = ConvexHull(X0, ϕ * X0 + δ * inputs) + if U isa ConstantInput + Ud = map(ui -> δ*ui ⊕ Eψ0, U) + elseif U isa VaryingInput + Ud = Vector{LazySet}(undef, length(Uset)) # TODO: use concrete type for Uset? how? + for (k, Uk) in enumerate(U) + Eψk = sih(Phi2Aabs * sih(A * Uk)) + Ud[k] = δ * Ui ⊕ Eψk end + Ud = VaryingInput(Ud) else - EPsi = sih(Phi2Aabs * sih(A * inputs)) - discretized_U = δ * inputs + EPsi - if approximation == "forward" - EOmegaPlus = sih(Phi2Aabs * sih((A * A) * X0)) - Ω0 = ConvexHull(X0, ϕ * X0 + discretized_U + EOmegaPlus) - elseif approximation == "backward" - EOmegaMinus = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) - Ω0 = ConvexHull(X0, ϕ * X0 + discretized_U + EOmegaMinus) - end + throw(ArgumentError("input of type $(typeof(U)) is not allwed")) end - if U isa ConstantInput - return DiscreteSystem(ϕ, Ω0, discretized_U) - else - discretized_U = [δ * Ui + sih(Phi2Aabs * sih(A * Ui)) for Ui in U] - return DiscreteSystem(ϕ, Ω0, discretized_U) - end + return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) end From 17f4dd911a59137bf1ddb7aea5e3846d0e92a85e Mon Sep 17 00:00:00 2001 From: mforets Date: Mon, 4 Mar 2019 18:40:37 -0300 Subject: [PATCH 04/24] update tests in systems and in discretization --- src/ReachSets/discretize.jl | 26 ++++---- test/Properties/unit_tune.jl | 4 +- test/ReachSets/unit_discretization.jl | 75 +++++++++++++---------- test/Reachability/models/bouncing_ball.jl | 2 +- test/Reachability/models/thermostat.jl | 4 +- test/Reachability/solve_continuous.jl | 30 ++++----- test/Systems/unit_System.jl | 20 +++--- test/runtests.jl | 8 ++- 8 files changed, 92 insertions(+), 77 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 58a5f987..20c46af0 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -1,7 +1,7 @@ const LDS = LinearDiscreteSystem const CLCDS = ConstrainedLinearControlDiscreteSystem -@inline I(T, n) = Matrix{eltype(A)}(I, n, n) +@inline Id(n) = Matrix{Float64}(I, n, n) """ discretize(𝑆, δ; [approximation], [exp_method], [sih_method]) @@ -399,7 +399,7 @@ function _discretize_firstorder(𝑆::InitialValueProblem, # compute exp(A*δ) ϕ = exp_Aδ(A, δ, exp_method=exp_method) - if islinear(𝑆) # inputdim(𝑆) == 0 + if inputdim(𝑆) == 0 α = κ * RX0 □α = Ballp(p, zeros(n), α) Ω0 = ConvexHull(X0, ϕ * X0 ⊕ □α) @@ -421,8 +421,8 @@ function _discretize_firstorder(𝑆::InitialValueProblem, # transformation of the inputs □β = Ballp(p, zeros(n), β) - Ud = map(u -> δ*u ⊕ □β, U) - return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) + Ud = map(u -> δ*u ⊕ □β, Uset) + return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) elseif Uset isa VaryingInput Ud = Vector{LazySet}(undef, length(Uset)) # TODO: use concrete type of Uset? how? @@ -444,7 +444,7 @@ function _discretize_firstorder(𝑆::InitialValueProblem, Ud[i] = δ*Ui ⊕ □β end Ud = VaryingInput(Ud) - return IVP(CLCDS(ϕ, I(typeof(ϕ), n), nothing, Ud), Ω0) + return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) end else throw(ArgumentError("this function only applies to linear or affine systems")) @@ -520,16 +520,16 @@ function _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousS Ω0 = copy(X0) # early return for homogeneous systems - if islinear(𝑆) + if inputdim(𝑆) == 0 return IVP(LDS(ϕ), Ω0) end # compute matrix to transform the inputs - Phi1Adelta = Φ₁_Aδ(A, δ, exp_method=exp_method) + Phi1Adelta = Φ₁(A, δ, exp_method=exp_method) U = inputset(𝑆) Ud = map(ui -> Phi1Adelta * ui, U) - return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) + return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0) end """ @@ -622,7 +622,7 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou ϕ = exp_Aδ(A, δ, exp_method=exp_method) # compute the transformation matrix to bloat the initial states - Phi2Aabs = Φ₂_Aδ(abs.(A), δ, exp_method=exp_method) + Phi2Aabs = Φ₂(abs.(A), δ, exp_method=exp_method) if approximation == "forward" Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use Eplus @@ -633,7 +633,7 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou end # early return for homogeneous systems - if islinear(𝑆) + if inputdim(𝑆) == 0 Ω0 = ConvexHull(X0, ϕ * X0 ⊕ Einit) return IVP(LDS(ϕ), Ω0) end @@ -647,15 +647,15 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou if U isa ConstantInput Ud = map(ui -> δ*ui ⊕ Eψ0, U) elseif U isa VaryingInput - Ud = Vector{LazySet}(undef, length(Uset)) # TODO: use concrete type for Uset? how? + Ud = Vector{LazySet}(undef, length(U)) # TODO: use concrete type for Uset? how? for (k, Uk) in enumerate(U) Eψk = sih(Phi2Aabs * sih(A * Uk)) - Ud[k] = δ * Ui ⊕ Eψk + Ud[k] = δ * Uk ⊕ Eψk end Ud = VaryingInput(Ud) else throw(ArgumentError("input of type $(typeof(U)) is not allwed")) end - return IVP(CLCDS(ϕ, I(typeof(A), n), nothing, Ud), Ω0) + return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0) end 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..4bf65c0a 100644 --- a/test/ReachSets/unit_discretization.jl +++ b/test/ReachSets/unit_discretization.jl @@ -4,92 +4,99 @@ 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, δ, approximation="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, δ, approximation="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, δ, approximation="firstorder") + @test inputdim(discr_sys_homog) == 0 + + # bloating, forward interpolation + discr_sys_homog = discretize(cont_sys_homog, δ, approximation="forward") + @test inputdim(discr_sys_homog) == 0 + + # bloating, backward interpolation + discr_sys_homog = discretize(cont_sys_homog, δ, approximation="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, δ, approximation="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, LazySets.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, δ, approximation="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, δ, approximation="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, δ, approximation="nobloating") Ui_d = inputset(discr_sys) @test length(Ui_d) == 3 for (i, inputs) in enumerate(discr_sys.s.U) if i == 1 @test dim(inputs) == 4 - @test isa(inputs, LinearMap) + @test isa(inputs, LazySets.LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*ones(4) && inputs.X.radius == 0.2 elseif i == 2 @test dim(inputs) == 4 - @test isa(inputs, LinearMap) + @test isa(inputs, LazySets.LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*2*ones(4) && inputs.X.radius == 0.2*2 else @test dim(inputs) == 4 - @test isa(inputs, LinearMap) + @test isa(inputs, LazySets.LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*3*ones(4) && inputs.X.radius == 0.2*3 end end # no bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, approx_model="nobloating", pade_expm=true) + discr_sys = discretize(cont_sys, δ, approximation="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..819f391e 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, B*U) # modes modes = [m1] diff --git a/test/Reachability/models/thermostat.jl b/test/Reachability/models/thermostat.jl index 9d1d34a8..a078d4d0 100644 --- a/test/Reachability/models/thermostat.jl +++ b/test/Reachability/models/thermostat.jl @@ -19,7 +19,7 @@ function thermostat() 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) + A, Matrix(1.0I, size(B, 1), size(B, 1)), inv, B*U) # mode off A = hcat(-c_a) @@ -27,7 +27,7 @@ function thermostat() 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) + A, Matrix(1.0I, size(B, 1), size(B, 1)), inv, B*U) # modes modes = [m_on, m_off] diff --git a/test/Reachability/solve_continuous.jl b/test/Reachability/solve_continuous.jl index a429ccae..4491f3dc 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,32 +63,32 @@ 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)) -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)) -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)) -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)) # 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, @@ -98,7 +98,7 @@ s = solve(ContinuousSystem(sparse(A), X0), # 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") From 1b238881f231582479146d7d82c887779764f8fd Mon Sep 17 00:00:00 2001 From: mforets Date: Mon, 4 Mar 2019 18:55:31 -0300 Subject: [PATCH 05/24] update BFFPSV18 --- .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 20 ++++--------------- .../ContinuousPost/BFFPSV18/check_property.jl | 13 +++--------- .../ContinuousPost/BFFPSV18/reach.jl | 13 +++--------- src/ReachSets/discretize.jl | 2 +- src/Utils/Utils.jl | 4 ++-- 5 files changed, 13 insertions(+), 39 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index ac7a9560..1815e983 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -24,7 +24,7 @@ end function options_BFFPSV18() return OptionSpec[ # general options - OptionSpec(:approx_model, "forward", domain=String, domain_check=( + OptionSpec(:approximation, "forward", domain=String, domain_check=( v -> v in ["forward", "backward", "firstorder", "nobloating"]), info="model for bloating/continuous time analysis"), OptionSpec(:algorithm, "explicit", domain=String, domain_check=( @@ -41,16 +41,10 @@ function options_BFFPSV18() "containing its indices"), # discretization options - OptionSpec(:lazy_sih, false, domain=Bool, + OptionSpec(:sih_method, "concrete", domain=String, 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(: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 +202,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] diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index b0356d2b..a87f9a98 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 + Δ = @timing discretize(system, options[:δ], approximation=options[:approximation], + exp_method=options[:exp_method], + sih_method=options[:sih_method]) Δ = matrix_conversion_lazy_explicit(Δ, 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..989ec50b 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -227,16 +227,9 @@ 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 + Δ = @timing discretize(system, options[:δ], approximation=options[:approximation], + exp_method=options[:exp_method], + sih_method=options[:sih_method]) Δ = matrix_conversion_lazy_explicit(Δ, options) return reach(Δ, invariant, options) end diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 20c46af0..96b34e5c 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -1,7 +1,7 @@ const LDS = LinearDiscreteSystem const CLCDS = ConstrainedLinearControlDiscreteSystem -@inline Id(n) = Matrix{Float64}(I, n, n) +@inline Id(n) = Matrix(1.0I, n, n) """ discretize(𝑆, δ; [approximation], [exp_method], [sih_method]) diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 6d887a39..d97e215d 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -516,10 +516,10 @@ function matrix_conversion(Δ, options; A_passed=nothing) return Δ end -# convert SparseMatrixExp to eplicit matrix +# convert SparseMatrixExp to explicit matrix function matrix_conversion_lazy_explicit(Δ, options) A = Δ.s.A - if !options[:lazy_expm] && options[:lazy_expm_discretize] + if options[:exp_method] == "lazy" # TODO: check this line info("Making lazy matrix exponential explicit...") @timing begin n = options[:n] From 94e75968ee06c0e127d619b71479837bda99543b Mon Sep 17 00:00:00 2001 From: mforets Date: Mon, 4 Mar 2019 19:28:23 -0300 Subject: [PATCH 06/24] fix tests --- .../ContinuousPost/BFFPSV18/check_property.jl | 3 ++- src/ReachSets/ContinuousPost/BFFPSV18/reach.jl | 2 ++ src/Utils/Utils.jl | 8 +++++++- src/solve.jl | 2 +- test/Reachability/solve_continuous.jl | 14 +++++++------- 5 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index a87f9a98..9788f694 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -182,9 +182,10 @@ function check_property(S::IVP{<:AbstractContinuousSystem}, # Time discretization # =================== info("Time discretization...") - Δ = @timing discretize(system, options[:δ], approximation=options[:approximation], + Δ = @timing discretize(S, options[:δ], approximation=options[:approximation], exp_method=options[:exp_method], sih_method=options[:sih_method]) + # TODO: remove this conversion and the option assume_sparse ? Δ = matrix_conversion_lazy_explicit(Δ, options) return check_property(Δ, property, options) end diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl index 989ec50b..6ebaa4c6 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -230,6 +230,8 @@ function reach(system::IVP{<:AbstractContinuousSystem}, Δ = @timing discretize(system, options[:δ], approximation=options[:approximation], exp_method=options[:exp_method], sih_method=options[:sih_method]) + + # TODO: remove this conversion and the option assume_sparse ? Δ = matrix_conversion_lazy_explicit(Δ, options) return reach(Δ, invariant, options) end diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index d97e215d..e8ffcf6c 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -519,7 +519,10 @@ end # convert SparseMatrixExp to explicit matrix function matrix_conversion_lazy_explicit(Δ, options) A = Δ.s.A - if options[:exp_method] == "lazy" # TODO: check this line + + #= + old, TODO: remove or review if we want to have this option again + if !options[:lazy_expm] && options[:lazy_expm_discretize] info("Making lazy matrix exponential explicit...") @timing begin n = options[:n] @@ -535,6 +538,9 @@ function matrix_conversion_lazy_explicit(Δ, options) else B = nothing end + =# + + B = nothing return matrix_conversion(Δ, options; A_passed=B) end diff --git a/src/solve.jl b/src/solve.jl index d1730d6d..936e3f96 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(CLCCS(loc.A, Matrix(1.0A, size(loc.A)), nothing, loc.U), X0.X), options_copy, op=opC, invariant=source_invariant) diff --git a/test/Reachability/solve_continuous.jl b/test/Reachability/solve_continuous.jl index 4491f3dc..bddea11f 100644 --- a/test/Reachability/solve_continuous.jl +++ b/test/Reachability/solve_continuous.jl @@ -70,28 +70,28 @@ s = solve(IVP(LCS(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(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(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(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(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)) # =============================== From a13294ce401de86df2d422c7163619b46714cf95 Mon Sep 17 00:00:00 2001 From: mforets Date: Mon, 4 Mar 2019 21:50:50 -0300 Subject: [PATCH 07/24] remove usage of ContinuousSystem and DiscreteSystem and fix tests --- docs/src/lib/discretize.md | 13 ------------- src/Reachability.jl | 6 ++++++ src/Transformations/Transformations.jl | 7 +++++-- src/Utils/systems.jl | 16 +++++++++------- src/solve.jl | 2 +- 5 files changed, 21 insertions(+), 23 deletions(-) diff --git a/docs/src/lib/discretize.md b/docs/src/lib/discretize.md index a651ca3c..5c77b7bf 100644 --- a/docs/src/lib/discretize.md +++ b/docs/src/lib/discretize.md @@ -12,16 +12,3 @@ CurrentModule = Reachability.ReachSets ```@docs discretize ``` - -## Dense-time reachability - -```@docs -discr_bloat_firstorder -discr_bloat_interpolation -``` - -## Discrete-time reachability - -```@docs -discr_no_bloat -``` 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/systems.jl b/src/Utils/systems.jl index 38cefd72..72c7cfda 100644 --- a/src/Utils/systems.jl +++ b/src/Utils/systems.jl @@ -15,6 +15,7 @@ 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) @@ -145,7 +146,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 +155,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 +169,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 +193,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 936e3f96..848837e0 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!(IVP(CLCCS(loc.A, Matrix(1.0A, size(loc.A)), nothing, loc.U), X0.X), + reach_tube = solve!(IVP(CLCCS(loc.A, Matrix(1.0I, size(loc.A)), nothing, ConstantInput(loc.U)), X0.X), options_copy, op=opC, invariant=source_invariant) From 94244abe66e7bdd1dddd33a541d151c64989ab3f Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 10:40:35 -0300 Subject: [PATCH 08/24] use loc for issue 491 --- .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 2 +- src/Utils/Utils.jl | 25 +------------------ src/solve.jl | 2 +- test/Reachability/models/bouncing_ball.jl | 2 +- test/Reachability/models/thermostat.jl | 8 +++--- 5 files changed, 8 insertions(+), 31 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index 1815e983..8f46535c 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -42,7 +42,7 @@ function options_BFFPSV18() # discretization options OptionSpec(:sih_method, "concrete", domain=String, - info="use a lazy symmetric interval hull in discretization?"), + 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, diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index e8ffcf6c..a3c9d1d6 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -518,30 +518,7 @@ end # convert SparseMatrixExp to explicit matrix function matrix_conversion_lazy_explicit(Δ, options) - A = Δ.s.A - - #= - old, TODO: remove or review if we want to have this option again - 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 - =# - - B = nothing - return matrix_conversion(Δ, options; A_passed=B) + return matrix_conversion(Δ, options; A_passed=Δ.s.A) end end # module diff --git a/src/solve.jl b/src/solve.jl index 848837e0..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!(IVP(CLCCS(loc.A, Matrix(1.0I, size(loc.A)), nothing, ConstantInput(loc.U)), X0.X), + reach_tube = solve!(IVP(loc, X0.X), options_copy, op=opC, invariant=source_invariant) diff --git a/test/Reachability/models/bouncing_ball.jl b/test/Reachability/models/bouncing_ball.jl index 819f391e..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(1.0I, 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 a078d4d0..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(1.0I, 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(1.0I, size(B, 1), size(B, 1)), inv, B*U) + m_off = ConstrainedLinearControlContinuousSystem(A, I1, inv, ConstantInput(B*U)) # modes modes = [m_on, m_off] From efd3c6eb373661daad34ddaa2e004886264c2393 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 10:56:02 -0300 Subject: [PATCH 09/24] remove matrix conversion lazy explicit function --- src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl | 2 +- src/ReachSets/ContinuousPost/BFFPSV18/reach.jl | 2 +- src/Utils/Utils.jl | 9 ++------- 3 files changed, 4 insertions(+), 9 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index 9788f694..579106a4 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -186,6 +186,6 @@ function check_property(S::IVP{<:AbstractContinuousSystem}, exp_method=options[:exp_method], sih_method=options[:sih_method]) # TODO: remove this conversion and the option assume_sparse ? - Δ = matrix_conversion_lazy_explicit(Δ, options) + Δ = 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 6ebaa4c6..dfc93f93 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -232,7 +232,7 @@ function reach(system::IVP{<:AbstractContinuousSystem}, sih_method=options[:sih_method]) # TODO: remove this conversion and the option assume_sparse ? - Δ = matrix_conversion_lazy_explicit(Δ, options) + Δ = matrix_conversion(Δ, options) return reach(Δ, invariant, options) end diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index a3c9d1d6..99da711e 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -508,17 +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 explicit matrix -function matrix_conversion_lazy_explicit(Δ, options) - return matrix_conversion(Δ, options; A_passed=Δ.s.A) -end - end # module From c04fc23cd67e5fcc0ee2ef295be4b0f167b054da Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 11:09:42 -0300 Subject: [PATCH 10/24] use algorithm for discretize kwarg --- .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 2 +- src/ReachSets/discretize.jl | 36 +++++++++---------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index 8f46535c..b2c7f00a 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -24,7 +24,7 @@ end function options_BFFPSV18() return OptionSpec[ # general options - OptionSpec(:approximation, "forward", domain=String, domain_check=( + OptionSpec(:discretization, "forward", domain=String, domain_check=( v -> v in ["forward", "backward", "firstorder", "nobloating"]), info="model for bloating/continuous time analysis"), OptionSpec(:algorithm, "explicit", domain=String, domain_check=( diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 96b34e5c..0383fe13 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -4,7 +4,7 @@ const CLCDS = ConstrainedLinearControlDiscreteSystem @inline Id(n) = Matrix(1.0I, n, n) """ - discretize(𝑆, δ; [approximation], [exp_method], [sih_method]) + discretize(𝑆, δ; [algorithm], [exp_method], [sih_method]) Apply an approximation model to `S` obtaining a discrete initial value problem. @@ -13,8 +13,8 @@ Apply an approximation model to `S` obtaining a discrete initial value problem. - `𝑆` -- initial value problem for a continuous affine ODE with non-deterministic inputs - `δ` -- step size -- `approximation` -- the method to compute the approximation model for the - discretization, choose among: +- `algorithm` -- the algorithm used to compute the approximation model for + the discretization, choose among: - `"forward"` -- use forward-time interpolation - `"backward"` -- use backward-time interpolation @@ -61,7 +61,7 @@ together with the coefficient matrix ``ϕ = e^{Aδ}`` and a transformed set of inputs if `U` is non-empty. In the literature, the method to obtain `Ω₀` is called the *approximation model* -and different alternatives have been proposed. See the argument `approximation` +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. @@ -71,7 +71,7 @@ discretized system. 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 -`approximation="nobloating"` for this setting. +`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 @@ -81,19 +81,19 @@ over vectors when needed, `e^{δA} v` for each `v`. Use the option """ function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; - approximation::String="forward", + algorithm::String="forward", exp_method::String="base", sih_method::String="lazy") - if approximation in ["forward", "backward"] - return _discretize_interpolation(𝑆, δ, approximation=approximation, + if algorithm in ["forward", "backward"] + return _discretize_interpolation(𝑆, δ, algorithm=algorithm, exp_method=exp_method, sih_method=sih_method) - elseif approximation == "firstorder" + elseif algorithm == "firstorder" return _discretize_firstorder(𝑆, δ, exp_method=exp_method) - elseif approximation == "nobloating" + elseif algorithm == "nobloating" return _discretize_nobloating(𝑆, δ, exp_method=exp_method) else - throw(ArgumentError("the approximation model $approximation is unknown")) + throw(ArgumentError("the algorithm $algorithm is unknown")) end end @@ -533,7 +533,7 @@ function _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousS end """ - _discretize_interpolation(𝑆, δ; [approximation], [exp_method], [sih_method]) + _discretize_interpolation(𝑆, δ; [algorithm], [exp_method], [sih_method]) Compute bloating factors using forward or backward interpolation. @@ -541,8 +541,8 @@ Compute bloating factors using forward or backward interpolation. - `𝑆` -- a continuous system - `δ` -- step size -- `approximation` -- choose the approximation model among `"forward"` and - `"backward"` +- `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: @@ -603,7 +603,7 @@ Heidelberg, 2011. """ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; - approximation::String="forward", + algorithm::String="forward", exp_method::String="base", sih_method::String="lazy") @@ -624,12 +624,12 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou # compute the transformation matrix to bloat the initial states Phi2Aabs = Φ₂(abs.(A), δ, exp_method=exp_method) - if approximation == "forward" + if algorithm == "forward" Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use Eplus - elseif approximation == "backward" + elseif algorithm == "backward" Einit = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) # use Eminus else - throw(ArgumentError("the method $approximation is unknown")) + throw(ArgumentError("the algorithm $approximation is unknown")) end # early return for homogeneous systems From 4fccb7a35d65b62c750e4af85a8c2b216b28b5e1 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 11:18:47 -0300 Subject: [PATCH 11/24] change option's name to discretization --- .../ContinuousPost/BFFPSV18/check_property.jl | 6 +++--- .../ContinuousPost/BFFPSV18/reach.jl | 2 +- test/ReachSets/unit_discretization.jl | 20 +++++++++---------- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl index 579106a4..c2fb6962 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -182,9 +182,9 @@ function check_property(S::IVP{<:AbstractContinuousSystem}, # Time discretization # =================== info("Time discretization...") - Δ = @timing discretize(S, options[:δ], approximation=options[:approximation], - exp_method=options[:exp_method], - sih_method=options[:sih_method]) + Δ = @timing discretize(S, options[:δ], algorithm=options[:discretization], + exp_method=options[:exp_method], + sih_method=options[:sih_method]) # TODO: remove this conversion and the option assume_sparse ? Δ = matrix_conversion(Δ, options) return check_property(Δ, property, options) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl index dfc93f93..13d7c188 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -227,7 +227,7 @@ function reach(system::IVP{<:AbstractContinuousSystem}, # Time discretization # =================== info("Time discretization...") - Δ = @timing discretize(system, options[:δ], approximation=options[:approximation], + Δ = @timing discretize(system, options[:δ], algorithm=options[:discretization], exp_method=options[:exp_method], sih_method=options[:sih_method]) diff --git a/test/ReachSets/unit_discretization.jl b/test/ReachSets/unit_discretization.jl index 4bf65c0a..7bbadd13 100644 --- a/test/ReachSets/unit_discretization.jl +++ b/test/ReachSets/unit_discretization.jl @@ -8,11 +8,11 @@ let # preventing global scope δ = 0.01 # no bloating, use default options - discr_sys_homog = discretize(cont_sys_homog, δ, approximation="nobloating") + 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, δ, approximation="nobloating", exp_method="pade") + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="nobloating", exp_method="pade") @test inputdim(discr_sys_homog) == 0 # bloating, default options @@ -20,15 +20,15 @@ let # preventing global scope @test inputdim(discr_sys_homog) == 0 # bloating, first order - discr_sys_homog = discretize(cont_sys_homog, δ, approximation="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, δ, approximation="forward") + 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, δ, approximation="backward") + discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="backward") @test inputdim(discr_sys_homog) == 0 # =============================================================== @@ -39,7 +39,7 @@ let # preventing global scope cont_sys = IVP(CLCCS(A, B, nothing, U), X0) # no bloating - discr_sys = discretize(cont_sys, δ, approximation="nobloating") + discr_sys = discretize(cont_sys, δ, algorithm="nobloating") @test inputdim(discr_sys) == 4 inputs = next_set(inputset(discr_sys)) @@ -48,7 +48,7 @@ let # preventing global scope @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, δ, approximation="nobloating", exp_method="pade") + discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade") @test inputdim(discr_sys) == 4 # bloating, use scaling and squaring method @@ -63,7 +63,7 @@ let # preventing global scope discr_sys = discretize(cont_sys, δ, exp_method="pade") @test inputdim(discr_sys) == 4 - discr_sys = discretize(cont_sys, δ, approximation="firstorder") + discr_sys = discretize(cont_sys, δ, algorithm="firstorder") @test inputdim(discr_sys) == 4 # =================================================================== @@ -73,7 +73,7 @@ let # preventing global scope cont_sys = IVP(CLCCS(A, Matrix(1.0I, 4, 4), nothing, VaryingInput(Ui)), X0) # no bloating - discr_sys = discretize(cont_sys, δ, approximation="nobloating") + 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) @@ -93,7 +93,7 @@ let # preventing global scope end # no bloating, use Pade approximation - discr_sys = discretize(cont_sys, δ, approximation="nobloating", exp_method="pade") + discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade") # bloating discr_sys = discretize(cont_sys, δ, exp_method="base") From 275860e94e6f0e15e5e5e6eeacc3c0d25c903e0b Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 11:43:50 -0300 Subject: [PATCH 12/24] revision --- .../ContinuousPost/BFFPSV18/BFFPSV18.jl | 1 - .../ContinuousPost/BFFPSV18/check_property.jl | 1 - src/ReachSets/ContinuousPost/BFFPSV18/reach.jl | 1 - src/ReachSets/discretize.jl | 18 +++++++++--------- test/ReachSets/unit_discretization.jl | 8 ++++---- 5 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index b2c7f00a..445e7feb 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -325,7 +325,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 c2fb6962..2736a12b 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/check_property.jl @@ -185,7 +185,6 @@ function check_property(S::IVP{<:AbstractContinuousSystem}, Δ = @timing discretize(S, options[:δ], algorithm=options[:discretization], exp_method=options[:exp_method], sih_method=options[:sih_method]) - # TODO: remove this conversion and the option assume_sparse ? Δ = 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 13d7c188..9ee0e64c 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/reach.jl @@ -231,7 +231,6 @@ function reach(system::IVP{<:AbstractContinuousSystem}, exp_method=options[:exp_method], sih_method=options[:sih_method]) - # TODO: remove this conversion and the option assume_sparse ? Δ = matrix_conversion(Δ, options) return reach(Δ, invariant, options) end diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 0383fe13..13f4a216 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -33,7 +33,7 @@ Apply an approximation model to `S` obtaining a discrete initial value problem. the evaluation of the action of the matrix exponential using the `expmv` implementation from `Expokit` -- `sih_method` -- (optional, default: `"lazy"`) the method used to take the +- `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 @@ -83,7 +83,7 @@ function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; algorithm::String="forward", exp_method::String="base", - sih_method::String="lazy") + sih_method::String="concrete") if algorithm in ["forward", "backward"] return _discretize_interpolation(𝑆, δ, algorithm=algorithm, @@ -425,7 +425,7 @@ function _discretize_firstorder(𝑆::InitialValueProblem, return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) elseif Uset isa VaryingInput - Ud = Vector{LazySet}(undef, length(Uset)) # TODO: use concrete type of Uset? how? + Ud = Vector{LazySet}(undef, length(Uset)) for (i, Ui) in enumerate(Uset) RU = norm(Ui, p) @@ -555,7 +555,7 @@ Compute bloating factors using forward or backward interpolation. the evaluation of the action of the matrix exponential using the `expmv` implementation from `Expokit` -- `sih_method` -- (optional, default: `"lazy"`) the method used to take the +- `sih_method` -- (optional, default: `"concrete"`) the method used to take the symmetric interval hull operation, choose among: - `"concrete"` -- compute the full symmetric interval hull @@ -605,12 +605,12 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou δ::Float64; algorithm::String="forward", exp_method::String="base", - sih_method::String="lazy") + sih_method::String="concrete") - if sih_method == "lazy" - sih = SymmetricIntervalHull - elseif sih_method == "concrete" + if sih_method == "concrete" sih = symmetric_interval_hull + elseif sih_method == "lazy" + sih = SymmetricIntervalHull else throw(ArgumentError("the method $sih_method is unknown")) end @@ -647,7 +647,7 @@ function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuou if U isa ConstantInput Ud = map(ui -> δ*ui ⊕ Eψ0, U) elseif U isa VaryingInput - Ud = Vector{LazySet}(undef, length(U)) # TODO: use concrete type for Uset? how? + Ud = Vector{LazySet}(undef, length(U)) for (k, Uk) in enumerate(U) Eψk = sih(Phi2Aabs * sih(A * Uk)) Ud[k] = δ * Uk ⊕ Eψk diff --git a/test/ReachSets/unit_discretization.jl b/test/ReachSets/unit_discretization.jl index 7bbadd13..560bbae1 100644 --- a/test/ReachSets/unit_discretization.jl +++ b/test/ReachSets/unit_discretization.jl @@ -44,7 +44,7 @@ let # preventing global scope inputs = next_set(inputset(discr_sys)) @test dim(inputs) == 4 - @test isa(inputs, LazySets.LinearMap) + @test isa(inputs, LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == ones(4) && inputs.X.radius == 0.5 # no bloating, use Pade approximation @@ -79,15 +79,15 @@ let # preventing global scope for (i, inputs) in enumerate(discr_sys.s.U) if i == 1 @test dim(inputs) == 4 - @test isa(inputs, LazySets.LinearMap) + @test isa(inputs, LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*ones(4) && inputs.X.radius == 0.2 elseif i == 2 @test dim(inputs) == 4 - @test isa(inputs, LazySets.LinearMap) + @test isa(inputs, LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*2*ones(4) && inputs.X.radius == 0.2*2 else @test dim(inputs) == 4 - @test isa(inputs, LazySets.LinearMap) + @test isa(inputs, LinearMap) @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*3*ones(4) && inputs.X.radius == 0.2*3 end end From 4a0a19c1846cead306f70b50e1aa8a3aca1eeb7d Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 13:00:42 -0300 Subject: [PATCH 13/24] remove domain check when it is already used in discretize --- src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index 445e7feb..3344b7d9 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -24,11 +24,9 @@ end function options_BFFPSV18() return OptionSpec[ # general options - OptionSpec(:discretization, "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"), + OptionSpec(:algorithm, "explicit", domain=String, OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time], domain_check=(v -> v > 0.), info="time step"), OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=( From 1f40964ab61431ea3e45ba9270fd63586b7a20c7 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 13:28:04 -0300 Subject: [PATCH 14/24] remove leading underscore in discretizes algorithms --- docs/src/lib/discretize.md | 3 +++ src/ReachSets/discretize.jl | 20 ++++++++++---------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/docs/src/lib/discretize.md b/docs/src/lib/discretize.md index 5c77b7bf..8f03ca7f 100644 --- a/docs/src/lib/discretize.md +++ b/docs/src/lib/discretize.md @@ -11,4 +11,7 @@ CurrentModule = Reachability.ReachSets ```@docs discretize +discretize_interpolation +discretize_firstorder +discretize_nobloating ``` diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 13f4a216..9406d773 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -86,12 +86,12 @@ function discretize(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, sih_method::String="concrete") if algorithm in ["forward", "backward"] - return _discretize_interpolation(𝑆, δ, algorithm=algorithm, + return discretize_interpolation(𝑆, δ, algorithm=algorithm, exp_method=exp_method, sih_method=sih_method) elseif algorithm == "firstorder" - return _discretize_firstorder(𝑆, δ, exp_method=exp_method) + return discretize_firstorder(𝑆, δ, exp_method=exp_method) elseif algorithm == "nobloating" - return _discretize_nobloating(𝑆, δ, exp_method=exp_method) + return discretize_nobloating(𝑆, δ, exp_method=exp_method) else throw(ArgumentError("the algorithm $algorithm is unknown")) end @@ -296,7 +296,7 @@ function Φ₂(A, δ; exp_method="base") end """ - _discretize_firstorder(𝑆, δ; [p], [exp_method]) + discretize_firstorder(𝑆, δ; [p], [exp_method]) Apply a first-order approximation model to `S` obtaining a discrete initial value problem. @@ -375,13 +375,13 @@ 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 +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, +function discretize_firstorder(𝑆::InitialValueProblem, δ::Float64; p::Float64=Inf, exp_method::String="base") @@ -452,7 +452,7 @@ function _discretize_firstorder(𝑆::InitialValueProblem, end """ - _discretize_nobloating(𝑆, δ; [exp_method]) + discretize_nobloating(𝑆, δ; [exp_method]) Discretize a continuous system without bloating of the initial states, suitable for discrete-time reachability. @@ -506,7 +506,7 @@ The transformations are: Here we allow ``U`` to be a sequence of time varying non-deterministic input sets. """ -function _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, +function discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; exp_method::String="base") @@ -533,7 +533,7 @@ function _discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousS end """ - _discretize_interpolation(𝑆, δ; [algorithm], [exp_method], [sih_method]) + discretize_interpolation(𝑆, δ; [algorithm], [exp_method], [sih_method]) Compute bloating factors using forward or backward interpolation. @@ -601,7 +601,7 @@ uses ``E^-``. International Conference on Computer Aided Verification. Springer, Berlin, Heidelberg, 2011. """ -function _discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, +function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, δ::Float64; algorithm::String="forward", exp_method::String="base", From 960f89c20706c09fc16fd92ab9bbe2b756451e2e Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 13:44:43 -0300 Subject: [PATCH 15/24] comment deprecated use of systems --- src/Utils/systems.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Utils/systems.jl b/src/Utils/systems.jl index 72c7cfda..c5910317 100644 --- a/src/Utils/systems.jl +++ b/src/Utils/systems.jl @@ -15,6 +15,7 @@ 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) @@ -33,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] From 073fa47c1cb0960a6c032d22f92e39676574e6b2 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 14:10:05 -0300 Subject: [PATCH 16/24] fix OptionSpec for explicit --- src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl index 3344b7d9..94cc81ad 100644 --- a/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl +++ b/src/ReachSets/ContinuousPost/BFFPSV18/BFFPSV18.jl @@ -26,7 +26,8 @@ function options_BFFPSV18() # general options OptionSpec(:discretization, "forward", domain=String, info="model for bloating/continuous time analysis"), - OptionSpec(:algorithm, "explicit", domain=String, + OptionSpec(:algorithm, "explicit", domain=String, domain_check=( + v -> v in ["explicit", "wrap"]), info="algorithm backend"), OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time], domain_check=(v -> v > 0.), info="time step"), OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=( From cb9ab5b6530b6df0a913b9fbf843bab82d4e15c4 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 14:38:39 -0300 Subject: [PATCH 17/24] update docstrings --- src/ReachSets/discretize.jl | 35 +++++++++++++++-------------------- 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 9406d773..727d6a69 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -13,8 +13,9 @@ Apply an approximation model to `S` obtaining a discrete initial value problem. - `𝑆` -- initial value problem for a continuous affine ODE with non-deterministic inputs - `δ` -- step size -- `algorithm` -- the algorithm used to compute the approximation model for - the discretization, choose among: +- `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 @@ -134,6 +135,12 @@ function exp_Aδ(A::AbstractMatrix{Float64}, δ::Float64; exp_method="base") end end +@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]) @@ -190,21 +197,15 @@ Heidelberg, 2011. function Φ₁(A, δ; exp_method="base") n = size(A, 1) if exp_method == "base" - P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)])) + P = expmat(Matrix(Pmatrix(A, δ, n))) Φ₁_Aδ = P[1:n, (n+1):2*n] elseif exp_method == "lazy" - P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) + P = SparseMatrixExp(Pmatrix(A, δ, n)) Φ₁_Aδ = sparse(get_columns(P, (n+1):2*n)[1:n, :]) elseif exp_method == "pade" - P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) + P = padm(Pmatrix(A, δ, n)) Φ₁_Aδ = P[1:n, (n+1):2*n] else @@ -271,21 +272,15 @@ Heidelberg, 2011. function Φ₂(A, δ; exp_method="base") n = size(A, 1) if exp_method == "base" - P = expmat(Matrix([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)])) + P = expmat(Matrix(Pmatrix(A, δ, n))) Φ₂_Aδ = P[1:n, (2*n+1):3*n] elseif exp_method == "lazy" - P = SparseMatrixExp([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) + P = SparseMatrixExp(Pmatrix(A, δ, n)) Φ₂_Aδ = sparse(get_columns(P, (2*n+1):3*n)[1:n, :]) elseif exp_method == "pade" - P = padm([A*δ sparse(δ*I, n, n) spzeros(n, n); - spzeros(n, 2*n) sparse(δ*I, n, n); - spzeros(n, 3*n)]) + P = padm(Pmatrix(A, δ, n)) Φ₂_Aδ = P[1:n, (2*n+1):3*n] else From 1cf1fd17a8d1f197eb1dda8fef6a874e7884b132 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 15:45:51 -0300 Subject: [PATCH 18/24] update doctring --- src/ReachSets/discretize.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 727d6a69..4adbd037 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -63,8 +63,8 @@ set of inputs if `U` is non-empty. 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. +for available options. For the reference to the original papers, see the +docstring of each method `discretize_...`. 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 From bdcbf2572f7f1bd301bf62ad61ae2aa57b46ab90 Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Tue, 5 Mar 2019 17:15:37 -0300 Subject: [PATCH 19/24] Update src/ReachSets/discretize.jl Co-Authored-By: mforets --- src/ReachSets/discretize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 4adbd037..bc65dccc 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -318,7 +318,7 @@ value problem. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` +- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` is returned, - otherwise a constrained linear discrete systen is returned, `ConstrainedLinearControlDiscreteSystem`. From fcdb8aee500e6daf74c35b4b1208cf255b59087a Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Tue, 5 Mar 2019 17:18:27 -0300 Subject: [PATCH 20/24] Update src/ReachSets/discretize.jl Co-Authored-By: mforets --- src/ReachSets/discretize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index bc65dccc..9fbf0330 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -320,7 +320,7 @@ The initial value problem for a discrete system. In particular: - if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` is returned, -- otherwise a constrained linear discrete systen is returned, +- otherwise a constrained linear discrete system is returned, `ConstrainedLinearControlDiscreteSystem`. ### Algorithm From 3b77455db91e6fede2f0476574315b1fe81e0db7 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 17:18:48 -0300 Subject: [PATCH 21/24] simplify docs --- src/ReachSets/discretize.jl | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index 4adbd037..feda71d8 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -318,10 +318,8 @@ value problem. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` - is returned, -- otherwise a constrained linear discrete systen is returned, - `ConstrainedLinearControlDiscreteSystem`. +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ### Algorithm @@ -472,10 +470,8 @@ for discrete-time reachability. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` - is returned, -- otherwise a constrained linear discrete systen is returned, - `ConstrainedLinearControlDiscreteSystem`. +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ### Algorithm @@ -561,10 +557,8 @@ Compute bloating factors using forward or backward interpolation. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a linear discrete system, `LinearDiscreteSystem` - is returned, -- otherwise a constrained linear discrete systen is returned, - `ConstrainedLinearControlDiscreteSystem`. +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ## Algorithm From 7f331911e0b0366505e330b4208724ac496841ba Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 17:49:02 -0300 Subject: [PATCH 22/24] various revisions --- src/ReachSets/discretize.jl | 40 ++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index feda71d8..ef034bf7 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -318,7 +318,7 @@ value problem. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, - otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ### Algorithm @@ -334,8 +334,8 @@ 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_{V} = \\max_{u ∈ U} ‖u‖``. If only the support functions of ``\\mathcal{X}_0`` -and ``V`` are known, these values might be hard to compute for some norms. See +``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: @@ -346,7 +346,7 @@ where ``α = (e^{δ ‖A‖} - 1 - δ‖A‖)*(R_{\\mathcal{X}_0} + R_{U} / ‖A ``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, δ]``, that we denote ``R_{[0,δ]}(\\mathcal{X}_0)`` here, +interval ``[0, δ]``, which we denote ``R_{[0,δ]}(\\mathcal{X}_0)`` here, is included in ``Ω₀``: ```math @@ -362,6 +362,8 @@ d_H(Ω₀, R_{[0,δ]}(\\mathcal{X}_0)) ≤ \\frac{1}{4}(e^{δ ‖A‖} - 1) D_{\ 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 @@ -375,9 +377,9 @@ uses less conservative bounds. using support functions. Nonlinear Analysis: Hybrid Systems, 4(2), 250-262.* """ function discretize_firstorder(𝑆::InitialValueProblem, - δ::Float64; - p::Float64=Inf, - exp_method::String="base") + δ::Float64; + p::Float64=Inf, + exp_method::String="base") # unwrap coefficient matrix and initial states A, X0 = 𝑆.s.A, 𝑆.x0 @@ -438,6 +440,8 @@ function discretize_firstorder(𝑆::InitialValueProblem, end 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")) @@ -470,7 +474,7 @@ for discrete-time reachability. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, - otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ### Algorithm @@ -498,10 +502,10 @@ The transformations are: Here we allow ``U`` to be a sequence of time varying non-deterministic input sets. """ function discretize_nobloating(𝑆::InitialValueProblem{<:AbstractContinuousSystem}, - δ::Float64; - exp_method::String="base") + δ::Float64; + exp_method::String="base") - # unrwap coefficient matrix and initial states + # unwrap coefficient matrix and initial states A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) @@ -557,7 +561,7 @@ Compute bloating factors using forward or backward interpolation. The initial value problem for a discrete system. In particular: -- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, +- if the input system is homogeneous, a `LinearDiscreteSystem` is returned, - otherwise a `ConstrainedLinearControlDiscreteSystem` is returned. ## Algorithm @@ -571,10 +575,6 @@ with ``x(0) ∈ \\mathcal{X}_0``, ``u(t) ∈ U`` be the given continuous affine `𝑆`, where `U` is the set of non-deterministic inputs and ``\\mathcal{X}_0`` is the set of initial states. -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. - The transformations are: - ``Φ ← \\exp^{Aδ}``, @@ -604,7 +604,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous throw(ArgumentError("the method $sih_method is unknown")) end - # unrwap coefficient matrix and initial states + # unwrap coefficient matrix and initial states A, X0 = 𝑆.s.A, 𝑆.x0 # compute matrix ϕ = exp(Aδ) @@ -614,9 +614,9 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous Phi2Aabs = Φ₂(abs.(A), δ, exp_method=exp_method) if algorithm == "forward" - Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use Eplus + Einit = sih(Phi2Aabs * sih((A * A) * X0)) # use E⁺ elseif algorithm == "backward" - Einit = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) # use Eminus + Einit = sih(Phi2Aabs * sih((A * A * ϕ) * X0)) # use E⁻ else throw(ArgumentError("the algorithm $approximation is unknown")) end @@ -643,7 +643,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous end Ud = VaryingInput(Ud) else - throw(ArgumentError("input of type $(typeof(U)) is not allwed")) + throw(ArgumentError("input of type $(typeof(U)) is not allowed")) end return IVP(CLCDS(ϕ, Id(size(A, 1)), nothing, Ud), Ω0) From 406d38975e5ff2a7b5e4cd891e759fa3b485656e Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 17:54:32 -0300 Subject: [PATCH 23/24] use ContantInput constructor in discretize_firstorder --- src/ReachSets/discretize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index ef034bf7..a13d66f9 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -416,7 +416,7 @@ function discretize_firstorder(𝑆::InitialValueProblem, # transformation of the inputs □β = Ballp(p, zeros(n), β) - Ud = map(u -> δ*u ⊕ □β, Uset) + Ud = ConstantInput(δ*U ⊕ □β) return IVP(CLCDS(ϕ, Id(n), nothing, Ud), Ω0) elseif Uset isa VaryingInput From 4bd092022a6954bf8f48be47236453126c2d4be3 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 5 Mar 2019 17:56:00 -0300 Subject: [PATCH 24/24] use ConstantInput constructor in discretize_interpolation --- src/ReachSets/discretize.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/ReachSets/discretize.jl b/src/ReachSets/discretize.jl index a13d66f9..dbdc0256 100644 --- a/src/ReachSets/discretize.jl +++ b/src/ReachSets/discretize.jl @@ -634,7 +634,8 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous Ω0 = ConvexHull(X0, ϕ * X0 ⊕ δ*U0 ⊕ Eψ0 ⊕ Einit) if U isa ConstantInput - Ud = map(ui -> δ*ui ⊕ Eψ0, U) + Ud = ConstantInput(δ*U0 ⊕ Eψ0) + elseif U isa VaryingInput Ud = Vector{LazySet}(undef, length(U)) for (k, Uk) in enumerate(U) @@ -642,6 +643,7 @@ function discretize_interpolation(𝑆::InitialValueProblem{<:AbstractContinuous Ud[k] = δ * Uk ⊕ Eψk end Ud = VaryingInput(Ud) + else throw(ArgumentError("input of type $(typeof(U)) is not allowed")) end