Discretization

Interface for conservative time discretization methods.

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.

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.

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.

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.

See also Bogomolov et al. [BFF+18], Eq. (14).

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 [FLD+11]. The Backward method uses $E^-$.

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 [FLD+11]. The Forward method uses $E^+$.

FirstOrder{EM} <: AbstractApproximationModel

First-order approximation model.

Fields

• exp – (optional, default: BaseExp) exponentiation method

Algorithm

The transformations are [LG10]:

• $Φ ← \exp(Aδ)$,
• $Ω_0 ← CH(\mathcal{X}_0, Φ\mathcal{X}_0 ⊕ δU ⊕ B_ε)$

where $B_ε$ is the input ball of radius $ε$ centered in the origin.

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.

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.

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.

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.

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., Forets and Schilling [FS22] and references therein.

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

Φ₁(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 Frehse et al. [FLD+11].

Φ₂(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 Frehse et al. [FLD+11].

AbstractExpAlg

Abstract supertype for all exponentiation algorithms.

BaseExpAlg <: AbstractExpAlg

Matrix exponential using the scaling-and-squaring algorithm implemented in Julia Base.

Notes

The alias for this algorithm is :base.

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

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.

PadeExpAlg <: AbstractExpAlg

Matrix exponential for sparse matrices using Pade approximants.

Notes

The alias for this algorithm is :pade.

This algorithm requires to load the package Expokit.jl.