LGG09{N, AM, TN<:AbstractDirections} <: AbstractContinuousPost

Implementation of Girard - Le Guernic algorithm for reachability analysis using support functions.


  • δ – step-size of the discretization
  • approx_model – (optional, default: Forward) approximation model; see Notes below for possible options
  • template – (alias: dirs) struct that holds the directions (either lazily or concretely) for each support function evaluation defining the template
  • static – (optional, default: false) if true, use statically sized arrays
  • threaded – (optional, default: false) if true, use multi-threading to compute different template directions in parallel
  • sparse – (optional, default: false) if true, convert the matrix exponential obtained after discretization to a sparse matrix
  • cache – (optional, default: true) if true, use a cache for intermediate computations in the set recurrence loop
  • vars – (optional, default: missing) used to specify the variables instead of passing the template


The type fields are:

  • N – number type of the step-size
  • AM – type of the approximation model
  • TN – type of the abstract directions that define the template

The flag threaded=true specifies that the support functions along different directions should be computed in parallel. See the section on Multi-threading for details on how to setup the number of threads.


This is an implementation of the algorithm from [LGG09].

These methods are described at length in the dissertation [LG09].



In the following subsections we outline the method of [LGG09] to solve linear set-based recurrences using support functions, first the homogeneous case and then the inhomogeneous case without wrapping effect. We also discuss the special case of real eigenvalues.

Homogeneous case

Consider the set-based recurrence

\[X_{k+1} = \Phi X_k,\qquad 0 \leq k \leq N\]

where $\Phi \in \mathbb{R}^{n\times n}$ and $X_0 \subseteq \mathbb{R}^n$ are given. By unrwapping the recurrence, $X_k = \Phi^k X_0$ for all $k \geq 0$. Let $d \in \mathbb{R}^n$ be a given template direction. Using the property of support functions $\rho(d, A X) = \rho(A^T d, X)$ for any matrix $A$ and set $X$, we have that

\[ρ(d, X_k) = \rho(d, \Phi^k X_0) = \rho((\Phi^T)^k d, X_0).\]

In this way we are able to reason with the sequence $\{X_0, X_1, X_2, \ldots, X_N\}$ by evaluating the support function of the initial set $X_0$ along the directions $\{d, \Phi^T d, (\Phi^T)^2 d, \ldots, (\Phi^T)^N d\}$.

Inhomogeneous case

The inhomogeneous case generalizes the previous case by taking, at each step, the Minkowski sum with an element from the sequence $\{V_0, V_1, V_2, \ldots, V_N\}$:

\[X_{k+1} = \Phi X_k \oplus V_k,\qquad 0 \leq k \leq N.\]

Let us write such recurrence in the unrapped form,

\[\begin{aligned} \quad X_1 &= \Phi X_0 \oplus V_0 \\[1mm] \quad X_2 &= \Phi X_1 \oplus V_1 = \Phi^2 X_0 \oplus \Phi V_0 \oplus V_1 \\[1mm] \quad X_3 &= \Phi X_2 \oplus V_2 = \Phi^3 X_0 \oplus \Phi^2 V_0 \oplus \Phi V_1 \oplus V_2 \\[1mm] \quad &\vdots \\[1mm] \quad X_k &= \Phi^k X_0 \oplus \left( \bigoplus_{i=0}^{k-1} \Phi^{k-i-1} V_i \right) \end{aligned}\]

where the big Minkowski sum is just an abbreviation for $\Phi^{k-1} V_0 \oplus \Phi^{k-2} V_1 \oplus \Phi^{k-3} V_2 \oplus \ldots \oplus \Phi V_{k-2} \oplus V_{k-1}$.

Let $d \in \mathbb{R}^n$ be a given template direction. Using the additive property of support functions, $\rho(d, X \oplus Y) = \rho(d, X) + \rho(d, Y)$ for any sets $X$ and $Y$, we have that

\[\begin{aligned} \quad \rho(d, X_1) &= \rho(\Phi^T d, X_0) + \rho(d, V_0) \\[1mm] \quad \rho(d, X_2) &= \rho((\Phi^T)^2 d, X_0) + \rho(\Phi^T d, V_0) + \rho(d, V_1) \\[1mm] \quad \rho(d, X_3) &= \rho((\Phi^T)^3 d, X_0) + \rho((\Phi^T)^2 d, V_0) + \rho(\Phi^T d, V_1) + \rho(d, V_2) \\[1mm] \quad &\vdots \\[1mm] \quad \rho(d, X_k) &= \rho((\Phi^T)^k d, X_0) + \sum_{i=0}^{k-1} \rho( (\Phi^T)^{k-i-1} d, V_i). \end{aligned}\]

In a similar fashion to the homogeneous case, the method allows to efficiently reason about the the sequence $\{X_0, X_1, X_2, \ldots, X_N\}$ by evaluating the support function of the initial set $X_0$ and the input sets $\{V_k\}_k$ along the directions $\{d, \Phi^T d, (\Phi^T)^2 d, \ldots, (\Phi^T)^N d\}$. Implementation-wise, we update two sequences, one that accounts for the homogeneous term, and another sequence that accounts for the effect of the accumulated inputs.

Implementation details

The reach-set representation used is a TemplateReachSet, which stores the directions used (vector of vectors) and the support function evaluated at each direction (matrix, see below). The set representation, set(R::TemplateReachSet), is either a polyhedron in constraint form (HPolyhedron), or a polytope (HPolytope) if the directions are bounding, i.e. the template directions define a bounded set.

The computed support function values can accessed directly through the field sf::SN of each template reach-set. Here sf is an array view of type ::Matrix{N}(undef, length(dirs), NSTEPS): each row corresponds to one of the template directions and each column corresponds to a fixed iteration index $k \geq 0$.

If you use directions from the canonical basis of $\mathbb{R}^n$, it is recommended to define LazySets.Arrays.SingleEntryVector or "one-hot" arrays as they are commonly called, because there are methods that dispatch on such type of arrays efficiently.

Parallel formulation

The support functions of the sequence $\{X_k\}_k$ along different directions can be computed in parallel. Indeed, if $d_1$ and $d_2$ are two given template directions, two different processes can independently compute $\rho(d_1, X_k)$ and $\rho(d_2, X_k)$ for all $k = 0, 1, \ldots, N$ using the methods described above. Once both computations have finished, we can store the resulting support functions in the same array. Use the flag threaded=true to use this method.

Implementation-wise the function _reach_homog_dir_LGG09! spawns differen threads which populate the matrix ρℓ::Matrix{N}(undef, length(dirs), NSTEPS) with the computed values. Hence each thread computes a subset of distinct rows of ρℓ.

Real eigenvalues

If the spectrum of the state matrix only has real eigenvalues, the sequence of support functions can be computed efficiently if we work with a template consisting of eigenvectors of $\Phi^T$. This idea is described in [LGG09] and we recall it here for convenience.

The method stems from the fact that if $(\lambda, d)$ is an eigenvalue-eigenvector pair of the matrix $\Phi^T$, with $\lambda \in \mathbb{R}$, then $\Phi^T d = \lambda d$, and if we apply $\Phi^T$ on both sides of this identity, we get $(\Phi^T)^2 d = \Phi^T (\Phi^T d) = \Phi^T(\lambda d) = \lambda^2 d$. In more generality, it holds that $(\Phi^T)^k d = \lambda^k d$ for all $k \ge 1$. Applying this relation to the support function recurrence described above, we get for the general inhomogeneous and possibly time-varying inputs case:

\[\rho(d, X_k) = \rho(\lambda^k d, X_0) + \sum_{i=0}^{k-1} \rho(\lambda^{k-i-1} d, V_i).\]

To further simplify this formula, we analyze different cases of $\lambda$. If $\lambda = 0$, then $\rho(d, X_k) = \rho(d, V_k)$ for all $k \geq 1$, so we focus on either $\lambda$ being positive or negative. To further simplify the computation of $\rho(d, X_k)$, we can use the property $\rho(\lambda d, X) = \lambda \rho(d, X)$ if $\lambda \geq 0$. We now consider the cases $\lambda > 0$ and $\lambda < 0$.

Case $\lambda > 0$. Then $\lambda^k > 0$ for all $k \geq 1$, and

\[\rho(d, X_k) = \lambda^k \rho(d, X_0) + \sum_{i=0}^{k-1} \lambda^{k-i-1} \rho(d, V_i).\]

We are left with evaluating the support function only at $\rho(d, X_0)$ and $\rho(d, V_i)$ to construct the full sequence $\{\rho(d, X_k)\}_{k}$. Moreover, if the $V_i$'s are constant we can extract them from the right-hand side sum and use that

\[\sum_{i=0}^{k-1} \lambda^{k-i-1} = 1 + \lambda + \ldots + \lambda^{k-1} = \dfrac{1 - \lambda^k}{1 - \lambda}.\]

Case $\lambda < 0$. Since $\lambda^k = (-1)^k (-\lambda)^k$ and $\lambda < 0$, then $\lambda^k$ is positive if $k$ is even, otherwise it is negative. So we can write:

\[\rho(d, X_k) = (-\lambda)^k \rho((-1)^k d, X_0) + \sum_{i=0}^{k-1} (-\lambda)^{k-i-1} \rho((-1)^{k-i-1} d, V_i).\]

The main difference between this case and the previus one is that now we have to evaluate support functions $\rho(\pm d, X_0)$ and $\rho(\pm d, V_i)$. Again, simplification takes place if the $V_i$'s are constant and such special case is considered in the implementation.