Discretization
Discretize API
ReachabilityAnalysis.DiscretizationModule
— ModuleInterface for conservative time discretization methods.
LinearAlgebra.normalize
— Functionnormalize(system::AbstractSystem)
Transform a mathematical system to a normalized (or canonical) form.
Input
system
– system; it can be discrete or continuous
Output
Either the same system if it already conforms to a canonical form, or a new system otherwise.
Notes
The normalization procedure consists of transforming a given system type into a "canonical" format that is used internally. More details are given below.
Algorithm
The implementation of normalize
exploits MathematicalSystems
's' types, which carry information about the problem as a type parameter.
Homogeneous ODEs of the form $x' = Ax, x ∈ \mathcal{X}$ are canonical if the associated problem is a ConstrainedLinearContinuousSystem
and A
is a matrix. This type does not handle non-deterministic inputs.
Note that a LinearContinuousSystem
does not consider constraints on the state-space (such as an invariant); to specify state constraints, use a ConstrainedLinearContinuousSystem
. If the passed system is a LinearContinuousSystem
(i.e. no constraints) then the normalization fixes a universal set (Universe
) as the constraint set.
The generalization to canonical systems with constraints and possibly time-varying non-deterministic inputs is considered next. These systems are of the form $x' = Ax + u, u ∈ \mathcal{U}, x ∈ \mathcal{X}$. The system type is ConstrainedLinearControlContinuousSystem
, where A
is a matrix, X
is a set and U
is an input, that is, any concrete subtype of AbstractInput
.
If U
is not given as an input, normalization accepts either a LazySet
, or a vector of LazySet
s. In these cases, the sets are wrapped around an appropriate concrete input type.
If the system does not conform to a canonical form, the implementation tries to make the transformation; otherwise an error is thrown. In particular, ODEs of the form $x' = Ax + Bu$ are mapped into $x' = Ax + u, u ∈ B\mathcal{U}$, where now $u$ has the same dimensions as $x$.
The transformations described above are analogous in the discrete case, i.e. $x_{k+1} = A x_k$ and $x_{k+1} = Ax_{k} + u_k, u_k ∈ \mathcal{U}, x_k ∈ \mathcal{X}$ for the linear and affine cases respectively.
normalize(ivp::InitialValueProblem)
Transform an initial-value problem into a normalized (or canonical) form.
Input
ivp
– initial-value problem
Output
Either the same initial-value problem if it already conforms to a canonical form, or a new one otherwise.
Notes
This function extends normalize
for initial-value problems.
MathematicalSystems.discretize
— Functiondiscretize(ivp::IVP, δ, alg::AbstractApproximationModel)
Set-based conservative discretization of a continuous-time initial value problem into a discrete-time problem.
Input
ivp
– initial value problem for a linear ODE in canonical form (seeNotes
below)δ
– step sizealg
– algorithm used to compute the approximation model
Output
The initial value problem of a discrete system.
Notes
Different approximation algorithms and their respective options are described in the docstring of each method. Here is a list of all the available approximation models:
julia> subtypes(ReachabilityAnalysis.DiscretizationModule.AbstractApproximationModel)
9-element Vector{Any}:
Backward
CorrectionHull
FirstOrder
FirstOrderZonotope
Forward
ForwardBackward
NoBloating
SecondOrderddt
StepIntersect
Initial-value problems considered in this function are of the form
\[x' = Ax(t) + u(t),\qquad x(0) ∈ \mathcal{X}_0,\qquad (1)\]
and where $u(t) ∈ U(k)$ add where $\{U(k)\}_k$ is a sequence of sets of non-deterministic inputs and $\mathcal{X}_0$ is the set of initial states. Other problems, e.g. $x' = Ax(t) + Bu(t)$ can be brought to the canonical form with the function normalize
.
For references to the original papers introducing each algorithm, see the docstrings, e.g. ?Forward
.
Approximation models
ReachabilityAnalysis.NoBloatingModule.NoBloating
— TypeNoBloating{EM, SO, IT} <: AbstractApproximationModel
No bloating, or discrete-time, approximation model.
Fields
exp
– exponentiation methodsetops
– set operations methodinv
– (optional, default:false
) iftrue
, assume that the state matrix is invertible and use its inverse in theΦ
functions
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$
- $Ω_0 ← \mathcal{X}_0$
- $V(k) ← Φ₁(A, δ)U(k)$, $k ≥ 0$.
The function $Φ₁(A, δ)$ is defined in Φ₁
. We allow $U$ to be a sequence of time varying non-deterministic input sets.
See also Eqs.(14) in [BFFPSV18].
ReachabilityAnalysis.ForwardModule.Forward
— TypeForward{EM, SO, SI, IT, BT} <: AbstractApproximationModel
Forward approximation model.
Fields
exp
– exponentiation methodsetops
– set operations methodsih
– symmetric interval hullinv
– (optional, default:false
) iftrue
, assume that the state matrix is invertible and use its inverse in theΦ
functionsbackend
– (optional, default:nothing
) used if the algorithm needs to apply concrete polyhedral computations
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$,
- $Ω_0 ← CH(\mathcal{X}_0, Φ\mathcal{X}_0 ⊕ δU(0) ⊕ E_ψ(U(0), δ) ⊕ E^+(\mathcal{X}_0, δ))$,
- $V(k) ← δU(k) ⊕ E_ψ(U(k), δ)$.
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 [FRE11]. The Backward
method uses $E^-$.
ReachabilityAnalysis.BackwardModule.Backward
— TypeBackward{EM, SO, SI, IT, BT} <: AbstractApproximationModel
Backward approximation model.
Fields
exp
– exponentiation methodsetops
– set operations methodsih
– symmetric interval hullinv
– (optional, default:false
) iftrue
, assume that the state matrix is invertible and use its inverse in theΦ
functionsbackend
– (optional, default:nothing
) used if the algorithm needs to apply concrete polyhedral computations
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$,
- $Ω_0 ← CH(\mathcal{X}_0, Φ\mathcal{X}_0 ⊕ δU(0) ⊕ E_ψ(U(0), δ) ⊕ E^-(\mathcal{X}_0, δ))$,
- $V(k) ← δU(k) ⊕ E_ψ(U(k), δ)$.
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 [FRE11]. The Forward
method uses $E^+$.
ReachabilityAnalysis.FirstOrderModule.FirstOrder
— TypeFirstOrder{EM} <: AbstractApproximationModel
First-order approximation model.
Fields
exp
– (optional, default:BaseExp
) exponentiation method
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$,
- $Ω_0 ← CH(\mathcal{X}_0, Φ\mathcal{X}_0 ⊕ δU ⊕ B_ε)$
where $B_ε$ is the input ball of radius $ε$ centered in the origin.
Reference
C. Le Guernic and A. Girard: Reachability analysis of linear systems using support functions. NAHS 2010.
ReachabilityAnalysis.FirstOrderZonotopeModule.FirstOrderZonotope
— TypeFirstOrderZonotope{EM} <: AbstractApproximationModel
First order approximation model that works with zonotopes.
Fields
exp
– (optional, default:BaseExp
) exponentiation method
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$,
- $Ω_0 ← bloat(zono(CH(\mathcal{X}_0, Φ\mathcal{X}_0)), α + β)$
where $bloat(\mathcal{X}, ε)$ bloats the set $\mathcal{X}$ with the value $ε$, $zono(·)$ overapproximates its argument with a zonotope, and $α$ and $β$ are factors computed for the homogeneous system and the inputs, respectively.
Reference
A. Girard: Reachability of uncertain linear systems using zonotopes. HSCC 2005.
ReachabilityAnalysis.SecondOrderddtModule.SecondOrderddt
— TypeSecondOrderddt{EM, SO, SI, IT, BT} <: AbstractApproximationModel
Second-order approximation model used in the tool d/dt
. It can be used for overapproximation and underapproximation.
Fields
oa
– (optional, default:true
) flag to choose between overapproximation and underapproximationexp
– (optional, default:BaseExp
) exponentiation methodsetops
– (optional, default::lazy
) set operations methodsih
– (optional, default::concrete
) way to compute the symmetric interval hullinv
– (optional, default:false
) iftrue
, assume that the state matrix is invertible and use its inverse in theΦ
functionsbackend
– (optional, default:nothing
) used if the algorithm needs to apply concrete polyhedral computations
Algorithm
The transformations are:
- $Φ ← \exp(Aδ)$,
- $Ω_0 ← bloat(CH(\mathcal{X}_0, Φ\mathcal{X}_0), ε)$
where $bloat(\mathcal{X}, ε)$ bloats the set $\mathcal{X}$ with the value $ε$. If oa == false
, the bloating acts in an inverted way and shrinks the set.
Reference
E. Asarin, T. Dang, O. Maler, O. Bournez: Approximate reachability analysis of piecewise-linear dynamical systems. HSCC 2000.
ReachabilityAnalysis.CorrectionHullModule.CorrectionHull
— TypeCorrectionHull{EM} <: AbstractApproximationModel
Discretization using the correction hull of the matrix exponential.
Fields
exp
– exponentiation methodorder
– order of the Taylor series expansion of the matrix exponential
Algorithm
For the homogeneous case, this method implements the transformation:
\[Ω_0 = CH(X_0, e^{Aδ} X_0) ⊕ FX_0\]
where $F$ is the correction (interval) matrix.
For the inhomogeneous case, $x' = Ax + u$, $x ∈ X, u ∈ U$, implements $Ω_0 = CH(X_0, exp(Aδ) X0) ⊕ FX0$ where $F$ is the correction (interval) matrix.
In both cases, if $A$ is an interval matrix, the exponential is overapproximated using methods from IntervalMatrices.jl
.
ReachabilityAnalysis.StepIntersectModule.StepIntersect
— TypeStepIntersect{DM<:AbstractApproximationModel} <: AbstractApproximationModel
Approximation model that composes (intersecting) one step Forward of a given model with one step backward of the same model.
Fields
model
– approximation modelsetops
– set operations method
Notes
Let $x' = Ax$ with $x(0) ∈ X₀$. This methods consists of:
Compute the discretized system with step-size $δ$ obtaining $Ω0$ and the given approximation model
method
.Compute the (lazy) linear map $ΦX₀$. This set consists of the (exact) reachable states at the time point $δ$.
Apply the approximation model
method
with initial condition $ΦX₀$ one step backward in time, with the state transition matrix $-A$, call this reach-set $Ω₀₊$Intersect $Ω₀₋$ and $Ω₀₊$ and return such set. The intersection is done either lazily or concretely depending on the specified
setops
field.
Exponentiation
The state transition matrix of the linear ODE $x'(t) = Ax(t) + u(t)$ at time $δ > 0$ is $Φ = e^{Aδ}$, hence the algorithms usually require to compute exponential matrices. There are distinct ways to compute the matrix exponential $e^{Aδ}$ depending on the type of $A$ (see e.g. [HIH08]). The available methods can be used through the (unexported) function _exp
.
For high dimensional systems (typically n > 2000
), computing the matrix exponential is expensive hence it is preferable to compute the action of the matrix exponential over vectors when needed, that is, $e^{δA} v$ for each $v$. This method is particularly well-suited if A
is vert sparse. Use the option exp=:krylov
(or exp=:lazy
) for this purpose.
ReachabilityAnalysis.Exponentiation
— ModuleInterface to matrix exponential backends of different kinds.
Includes common integral computations arising in the discretization of linear differential equations using matrix methods. For applications see e.g. [1] and references therein.
[1] Conservative Time Discretization: A Comparative Study. Marcelo Forets and Christian Schilling (2022). Proceedings of the 17th International Conference on integrated Formal Methods (iFM), LNCS, vol. 13274, pp. 149-167. doi: 10.1007/978-3-031-07727-2_9, arXiv: 2111.01454.
ReachabilityAnalysis.Exponentiation._exp
— Function_exp(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp)
Compute the matrix exponential $e^{Aδ}$.
Input
A
– matrixδ
– step sizealg
– (optional, default:BaseExp
) the algorithm used to take the matrix exponential ofAδ
, possible options areBaseExp
,LazyExp
,PadeExp
andIntervalExp
see details in the Algorithm section below
Output
A matrix or a lazy wrapper of the matrix exponential, depending on alg
.
Algorithm
BaseExp
– (alias::base
) use Higham's scaling and squaring method implemented in Julia standard library; see?exp
for details; ifA
is a static array, uses the implementation inStaticArrays.jl
LazyExp
– (alias::lazy
) return a lazy wrapper type around the matrix exponential using the implementationLazySets.SparseMatrixExp
PadeExp
– (alias:pade
) apply Pade approximant method to compute the matrix exponential of a sparse matrix (requiresExpokit.jl
)IntervalExp
– (alias:interval
,taylor
) apply the Taylor series expansion of the matrix exponential with an interval remainder; works ifA
is an interval matrix
Notes
If the algorithm LazyExp
is used, actions of the matrix exponential are evaluated with an external library such as ExponentialUtilities.jl
or Expokit.jl
.
ReachabilityAnalysis.Exponentiation.Φ₁
— FunctionΦ₁(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp, [isinv]::Bool=false, [Φ]=nothing)
Evaluate the series
\[Φ₁(A, δ) = ∑_{i=0}^∞ \dfrac{δ^{i+1}}{(i+1)!}A^i,\]
where $A$ is a square matrix of order $n$ and $δ > 0$.
Input
A
– coefficients matrixδ
– step sizealg
– (optional, default:BaseExp
) the method used to take the matrix exponential of the coefficient matrix; see the documentation of_exp
for available optionsisinv
– (optional, default:false
) iftrue
, assume that the coefficients matrix is invertible and compute $A^{-1}$Φ
– (optional, default:nothing
) optionally pass the matrix exponential $e^{Aδ}$
Output
A matrix.
Algorithm
If $A$ is invertible, $Φ₁$ can be computed as
\[Φ₁(A, δ) = A^{-1}(e^{Aδ} - I_n),\]
where $I_n$ is the identity matrix of order $n$.
In the general case, implemented in this function, it can be computed as submatrices of the block matrix
\[P_{2n} = \exp \begin{pmatrix} Aδ && δI_n \\ 0 && 0 \end{pmatrix}.\]
It can be shown that
\[\exp(P_{2n}) = \begin{pmatrix} Φ(A, δ) && Φ₁(A, δ) \\ 0 && δI_n \end{pmatrix}.\]
where $Φ(A, δ) = e^{Aδ}$. In particular, Φ₁(A, δ) = P[1:n, (n+1):2*n]
. This method can be found in [FRE11].
ReachabilityAnalysis.Exponentiation.Φ₂
— FunctionΦ₂(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp, [isinv]::Bool=false, [Φ]=nothing)
Evaluate the series
\[Φ₂(A, δ) = ∑_{i=0}^∞ \dfrac{δ^{i+2}}{(i+2)!}A^i,\]
where $A$ is a square matrix of order $n$ and $δ > 0$.
Input
A
– coefficients matrixδ
– step sizealg
– (optional, default:BaseExp
) the method used to take the matrix exponential of the coefficient matrix; see the documentation of_exp
for available optionsisinv
– (optional, default:false
) iftrue
, assume that the coefficients matrix is invertible and compute $A^{-1}$Φ
– (optional, default:nothing
) optionally pass the matrix exponential $e^{Aδ}$
Output
A matrix.
Algorithm
If $A$ is invertible, $Φ₂$ can be computed as
\[Φ₂(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
\[P_{3n} = \exp \begin{pmatrix} Aδ && δI_n && 0 \\ 0 && 0 && δI_n \\ 0 && 0 && 0 \end{pmatrix}.\]
It can be shown that
\[\exp(P_{3n}) = \begin{pmatrix} Φ(A, δ) && Φ₁(A, δ) && Φ₂(A, δ) \\ 0 && I_n && δI_n \\ 0 && 0 && I_n \end{pmatrix}.\]
where $Φ(A, δ) = e^{Aδ}$. In particular, Φ₂ = P_{3n}[1:n, (2*n+1):3*n]
. This method can be found in [FRE11].
ReachabilityAnalysis.Exponentiation.AbstractExpAlg
— TypeAbstractExpAlg
Abstract supertype for all exponentiation algorithms.
ReachabilityAnalysis.Exponentiation.BaseExpAlg
— TypeBaseExpAlg <: AbstractExpAlg
Matrix exponential using the scaling and squaring algorithm implemented in Julia Base.
Notes
If the array is static, the method implemented in StaticArrays.jl
is applied. This algorithm admits the alias :base
.
ReachabilityAnalysis.Exponentiation.LazyExpAlg
— TypeLazyExpAlg <: AbstractExpAlg
Lazy wrapper for the matrix exponential, operations defined in LazySets.jl
.
Fields
m
– (optional, default:30
) size of the Krylov subspacetol
– (optional, default:1e-10
) tolerance
Notes
This algorithm admits the alias :lazy
and also :krylov
(using default options).
ReachabilityAnalysis.Exponentiation.IntervalExpAlg
— TypeIntervalExpAlg <: AbstractExpAlg
Matrix exponential using an interval enclosure of the Taylor series remainder.
Fields
order
– order of the Taylor series expansion of the matrix exponential
Notes
This method allows to overapproximate $exp(Aδ)$ with an interval matrix. It also accepts an interval matrix $A$ given as input. This method admits the alias :interval
and :taylor
.
ReachabilityAnalysis.Exponentiation.PadeExpAlg
— TypePadeExpAlg <: AbstractExpAlg
Matrix exponential for sparse matrices using Pade approximants.
Notes
Requires Expokit.jl
. This algorithm admits the alias :pade
.
References
- HIH08Higham, Nicholas J. Functions of matrices: theory and computation. Society for Industrial and Applied Mathematics, 2008.