# Discretization

## Discretize API

`ReachabilityAnalysis.DiscretizationModule`

— ModuleInterface for conservative time discretization methods.

`LinearAlgebra.normalize`

— Function`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 `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`

— Function`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`

.

## Approximation models

`ReachabilityAnalysis.NoBloatingModule.NoBloating`

— Type`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 Eqs.(14) in [BFFPSV18].

`ReachabilityAnalysis.ForwardModule.Forward`

— Type`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^-$.

`ReachabilityAnalysis.BackwardModule.Backward`

— Type`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^+$.

`ReachabilityAnalysis.FirstOrderModule.FirstOrder`

— Type`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.

`ReachabilityAnalysis.FirstOrderZonotopeModule.FirstOrderZonotope`

— Type`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.

`ReachabilityAnalysis.SecondOrderddtModule.SecondOrderddt`

— Type`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.

`ReachabilityAnalysis.CorrectionHullModule.CorrectionHull`

— Type`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`

.

`ReachabilityAnalysis.StepIntersectModule.StepIntersect`

— Type`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.

## 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 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`

.

`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].

`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].

`ReachabilityAnalysis.Exponentiation.AbstractExpAlg`

— Type`AbstractExpAlg`

Abstract supertype for all exponentiation algorithms.

`ReachabilityAnalysis.Exponentiation.BaseExpAlg`

— Type`BaseExpAlg <: AbstractExpAlg`

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

**Notes**

The alias for this algorithm is `:base`

.

`ReachabilityAnalysis.Exponentiation.LazyExpAlg`

— Type`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`

).

`ReachabilityAnalysis.Exponentiation.IntervalExpAlg`

— Type`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.

`ReachabilityAnalysis.Exponentiation.PadeExpAlg`

— Type`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`

.

## References

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