Discretize

discretize(𝑆, δ; [algorithm], [exp_method], [sih_method], [set_operations])

Apply an approximation model to S obtaining a discrete initial value problem.

Input

• 𝑆 – initial value problem for a continuous affine ODE with non-deterministic inputs

• δ – step size

• 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
  • "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: "concrete") the method used to take the symmetric interval hull operation; choose among:

  • "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

• set_operations – (optional, default: "lazy") set operations used for the discretized initial states and transformed inputs; choose among:

  • "lazy" – use lazy convex hull for the initial states and lazy linear map for the inputs
  • "zonotope" – use concrete zonotope operations (linear map and Minkowski sum) and return zonotopes for both the initial states and the inputs of the discretized system

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 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, Ω₀, 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 Ω₀ 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 algorithm 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 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 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 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.

discretize_interpolation(𝑆, δ; [algorithm], [exp_method], [sih_method])

Compute bloating factors using forward or backward interpolation.

Input

• 𝑆 – a continuous system

• δ – step size

• 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:

  • "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: "concrete") 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

• set_operations – (optional, default: "lazy") set operations used for the discretized initial states and transformed inputs; choose among:

  • "lazy" – use lazy convex hull for the initial states and lazy linear map for the inputs
  • "zonotope" – use concrete zonotope operations (linear map and Minkowski sum) and return zonotopes for both the initial states and the inputs of the discretized system

Output

The initial value problem for a discrete system. In particular:

• if the input system is homogeneous, a LinearDiscreteSystem is returned,
• otherwise a ConstrainedLinearControlDiscreteSystem is returned.

Algorithm

Let us define some notation. Let

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 transformations are:

• $Φ ← \exp^{Aδ}$,
• $Ω₀ ← ConvexHull(\mathcal{X}_0, Φ\mathcal{X}_0 ⊕ δU(0) ⊕ Eψ(U(0), δ) ⊕ E^+(\mathcal{X}_0, δ))$,
• $V ← δ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 [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.

discretize_firstorder(𝑆, δ; [p], [exp_method])

Apply a first-order approximation model to S obtaining a discrete initial value problem.

Input

• 𝑆 – 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, 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. In particular:

• if the input system is homogeneous, a LinearDiscreteSystem is returned,
• otherwise a ConstrainedLinearControlDiscreteSystem is returned.

Algorithm

Let us define some notation. Let

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.

Define $R_{\mathcal{X}_0} = \max_{x ∈ \mathcal{X}_0} ‖x‖$, D_{\mathcal{X}_0} = \max_{x, y ∈ \mathcal{X}_0} ‖x-y‖andR_{U} = \max_{u ∈ U} ‖u‖. If only the support functions of\mathcal{X}_0andUare known, these values might be hard to compute for some norms. SeeNotes` below for details.

Let $Ω₀$ be the set defined as:

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, δ]$, which we denote $R_{[0,δ]}(\mathcal{X}_0)$ here, is included in $Ω₀$:

Moreover, if $d_H(A, B)$ denotes the Hausdorff distance between the sets $A$ and $B$ in $\mathbb{R}^n$, then

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

discretize_nobloating(𝑆, δ; [exp_method])

Discretize a continuous system without bloating of the initial states, suitable for discrete-time reachability.

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:

  • "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. In particular:

• if the input system is homogeneous, a LinearDiscreteSystem is returned,
• otherwise a ConstrainedLinearControlDiscreteSystem is returned.

Algorithm

Let us define some notation. Let

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:

• $Φ ← \exp^{Aδ}$
• $Ω₀ ← \mathcal{X}_0$
• $V ← Φ₁(A, δ)U(k)$, where $Φ₁(A, δ)$ is defined in Φ₁.

Here we allow $U$ to be a sequence of time varying non-deterministic input sets.

discretize_interval_matrix(𝑆::InitialValueProblem, δ::Float64, order::Int)

Discretize an affine continuous system whose dynamics are described by an interval matrix.

Input

• 𝑆 – a continuous system

• δ – step size

• order – order for Taylor-series approximation, ≥ 2

• set_operations – (optional, default: "zonotope") set operations used for the discretized initial states and transformed inputs; choose among:

  • "lazy" – use lazy convex hull for the initial states and lazy linear map for the inputs
  • "zonotope" – use concrete zonotope operations (linear map and Minkowski sum) and return zonotopes for both the initial states and the inputs of the discretized system

Output

The discretized system. In particular:

• if the input system is homogeneous, a LinearDiscreteSystem is returned,
• otherwise a ConstrainedLinearControlDiscreteSystem is returned.

Algorithm

See the first two lines of Algorithm 1 in [1].

[1] M. Althoff, O. Stursberg, M. Buss. Reachability Analysis of Linear Systems with Uncertain Parameters and Inputs. CDC 2007.

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, 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

A matrix.

Φ₁(A, δ; [exp_method])

Compute the series

where $A$ is a square matrix of order $n$ and $δ ∈ \mathbb{R}_{≥0}$.

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, 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

A matrix.

Algorithm

We use the method from [1]. If $A$ is invertible, $Φ₁$ can be computed as

In the general case, implemented in this function, it can be computed as submatrices of the block matrix

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.

Φ₂(A, δ; [exp_method])

Compute the series

where $A$ is a square matrix of order $n$ and $δ ∈ \mathbb{R}_{≥0}$.

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, 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

A matrix.

Algorithm

We use the method from [1]. If $A$ is invertible, $Φ₂$ can be computed as

In the general case, implemented in this function, it can be computed as submatrices of the block matrix

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.