# Discretization

## Discretize API

LinearAlgebra.normalizeFunction
normalize(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 LazySets. 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.

source
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.

source
MathematicalSystems.discretizeFunction
discretize(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 (see Notes below)
• δ – step size
• alg – 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.

source

## Approximation models

ReachabilityAnalysis.NoBloatingModule.NoBloatingType
NoBloating{EM, SO, IT} <: AbstractApproximationModel

No bloating, or discrete-time, approximation model.

Fields

• exp – exponentiation method
• setops – set operations method
• inv – (optional, default: false) if true, 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.

source
ReachabilityAnalysis.ForwardModule.ForwardType
Forward{EM, SO, SI, IT, BT} <: AbstractApproximationModel

Forward approximation model.

Fields

• exp – exponentiation method
• setops – set operations method
• sih – symmetric interval hull
• inv – (optional, default: false) if true, assume that the state matrix is invertible and use its inverse in the Φ functions
• backend – (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^-$.

source
ReachabilityAnalysis.BackwardModule.BackwardType
Backward{EM, SO, SI, IT, BT} <: AbstractApproximationModel

Backward approximation model.

Fields

• exp – exponentiation method
• setops – set operations method
• sih – symmetric interval hull
• inv – (optional, default: false) if true, assume that the state matrix is invertible and use its inverse in the Φ functions
• backend – (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^+$.

source
ReachabilityAnalysis.FirstOrderModule.FirstOrderType
FirstOrder{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.

source
ReachabilityAnalysis.FirstOrderZonotopeModule.FirstOrderZonotopeType
FirstOrderZonotope{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.

source
ReachabilityAnalysis.SecondOrderddtModule.SecondOrderddtType
SecondOrderddt{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 underapproximation
• exp – (optional, default: BaseExp) exponentiation method
• setops – (optional, default: :lazy) set operations method
• sih – (optional, default: :concrete) way to compute the symmetric interval hull
• inv – (optional, default: false) if true, assume that the state matrix is invertible and use its inverse in the Φ functions
• backend – (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.

source
ReachabilityAnalysis.CorrectionHullModule.CorrectionHullType
CorrectionHull{EM} <: AbstractApproximationModel

Discretization using the correction hull of the matrix exponential.

Fields

• exp – exponentiation method
• order – 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.

source
ReachabilityAnalysis.StepIntersectModule.StepIntersectType
StepIntersect{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 model
• setops – 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.

source

## 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.ExponentiationModule

Interface 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.

source
ReachabilityAnalysis.Exponentiation._expFunction
_exp(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp)

Compute the matrix exponential $e^{Aδ}$.

Input

• A – matrix
• δ – step size
• alg – (optional, default: BaseExp) the algorithm used to take the matrix exponential of Aδ, possible options are BaseExp, LazyExp, PadeExp and IntervalExp (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's standard library; see ?exp for details; if A is a static array, uses the implementation in StaticArrays.jl

• LazyExp – (alias: :lazy) return a lazy wrapper around the matrix exponential using the implementation LazySets.SparseMatrixExp

• PadeExp – (alias: pade) apply the Padé approximant method to compute the matrix exponential of a sparse matrix (requires Expokit.jl)

• IntervalExp – (alias: interval, taylor) apply the Taylor series expansion of the matrix exponential with an interval remainder; works if A 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.

source
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 size
• alg – (optional, default: BaseExp) the method used to take the matrix exponential of the coefficient matrix; see the documentation of _exp for available options
• isinv – (optional, default: false) if true, 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].

source
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 size
• alg – (optional, default: BaseExp) the method used to take the matrix exponential of the coefficient matrix; see the documentation of _exp for available options
• isinv – (optional, default: false) if true, 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].

source
ReachabilityAnalysis.Exponentiation.LazyExpAlgType
LazyExpAlg <: AbstractExpAlg

Matrix exponential computed in a lazy fashion.

Fields

• m – (optional, default: 30) size of the Krylov subspace
• tol – (optional, default: 1e-10) tolerance

Notes

The aliases for this algorithm are :lazy and :krylov.

The operations are defined in the package LazySets.jl (SparseMatrixExp).

source
ReachabilityAnalysis.Exponentiation.IntervalExpAlgType
IntervalExpAlg <: 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

The aliases for this algorithm are :interval and :taylor.

This algorithm allows to overapproximate $exp(Aδ)$ with an interval matrix. It also accepts an interval matrix $A$ as input.

source

## References

• HIH08Higham, Nicholas J. Functions of matrices: theory and computation. Society for Industrial and Applied Mathematics, 2008.