pymor.reductors.era

Module Contents

class pymor.reductors.era.ERAReductor(data, sampling_time, force_stability=True, feedthrough=None)

Bases: ERAReductorBase

Eigensystem Realization Algorithm reductor.

Constructs a (reduced) realization from a sequence of Markov parameters \(h_i\), for \(i\in\{1,\,\dots,\,2s-1\}\), \(s\in\mathbb{N}\), by a (reduced) orthogonal factorization of the Hankel matrix of Markov parameters

\[\begin{split}H = \begin{bmatrix} h_1 & h_2 & \dots & h_s \\ h_2 & h_3 & \dots & h_{s+1}\\ \vdots & \vdots & \ddots & \vdots\\ h_s & h_{s+1} & \dots & h_{2s-1} \end{bmatrix}=U\Sigma V^T\in\mathbb{R}^{ps\times ms},\end{split}\]

where \(r\leq\min\{ms,ps\}\) is the reduced order. See [Kun78].

In order for the identified model to be stable, the Markov parameters decay substantially within \(s\) samples. Stability is enforced automatically through zero-padding and can be deactivated by setting force_stability=False.

For a large number of inputs and/or outputs, the factorization of the Hankel matrix can be accelerated by tangentially projecting the Markov parameters to reduce the dimension of the Hankel matrix, i.e.

\[\hat{h}_i = W_L^T h_i W_R,\]

where \(n_L \leq p\) and \(n_R \leq m\) are the number of left and right tangential directions and \(W_L \in \mathbb{R}^{p \times n_L}\) an \(W_R \in \mathbb{R}^{m \times n_R}\) are the left and right projectors, respectively. See [KG16].

data

NumPy array that contains the first \(n\) Markov parameters of an LTI system. Has to be one- or three-dimensional with either:

data.shape == (n,)

for scalar-valued Markov parameters or:

data.shape == (n, p, m)

for matrix-valued Markov parameters of dimension \(p\times m\), where \(m\) is the number of inputs and \(p\) is the number of outputs of the system.

sampling_time

A number that denotes the sampling time of the system (in seconds).

force_stability

Whether the Markov parameters are zero-padded to double the length in order to enforce Kung’s stability assumption. See [Kun78]. Defaults to True.

feedthrough

(Optional) Operator or NumPy array of shape (p, m). The zeroth Markov parameter that defines the feedthrough of the realization. Defaults to None.

Methods

error_bounds	Compute the error bounds for all possible reduction orders.
reduce	Construct a minimal realization.

error_bounds(num_left=None, num_right=None)

Compute the error bounds for all possible reduction orders.

Without tangential projection of the Markov parameters, the \(\mathcal{L}_2\)-error of the Markov parameters \(\epsilon\) is bounded by

\[\epsilon = \lVert G-G_r\rVert_{\mathcal{H}_2}= \left(\sum_{i = 1}^{2 s - 1} \lVert C_r A_r^{i - 1} B_r - h_i \rVert_F^2\right)^{1/2} \leq \sigma_{r + 1}(\mathcal{H}) \sqrt{r + p + m},\]

where \((A_r,B_r,C_r)\) is the reduced realization of order \(r\), \(h_i\in\mathbb{R}^{p\times m}\) is the \(i\)-th Markov parameter and \(\sigma_{r+1}(\mathcal{H})\) is the first neglected singular value of the Hankel matrix of Markov parameters.

With tangential projection, the bound is given by

\[\epsilon = \lVert G-G_r\rVert_{\mathcal{H}_2}= \left(\sum_{i = 1}^{2 s - 1} \lVert C_r A_r^{i - 1} B_r - h_i \rVert_F^2\right)^{1/2} \leq 4 \left( \sum_{i = n_L + 1}^p \sigma_i^2(\Theta_L) + \sum_{i = n_R + 1}^m \sigma_i^2(\Theta_R) \right) + 2 \sigma_{r + 1}(\mathcal{H}) \sqrt{r + n_L + n_R},\]

where \(\Theta_L,\,\Theta_R\) is the matrix of horizontally or vertically stacked Markov parameters, respectively. See [KG16] (Thm. 3.4) for details and note that the present bound is squared due to a typographical error in [Kun78] that was reported in [PS26].

reduce(r=None, tol=None, num_left=None, num_right=None)

Construct a minimal realization.

Parameters:

• r – Order of the reduced model if tol is None, maximum order if tol is specified.

• tol – Tolerance for the error bound.

• num_left – Number of left (output) directions for tangential projection.

• num_right – Number of right (input) directions for tangential projection.

Returns:

rom – Reduced-order LTIModel.

class pymor.reductors.era.ERAReductorBase(data, sampling_time, force_stability=True, feedthrough=None)

Bases: pymor.core.cache.CacheableObject

Basic Eigensystem Realization Algorithm reductor.

This class implements ROM construction from a orthogonal factorization of the Hankel matrix, as well as tangential projections of inputs and outputs. The actual factorization and error bounds/estimators are implemented by the subclasses.

data

NumPy array that contains the first \(n\) Markov parameters of an LTI system. Has to be one- or three-dimensional with either:

data.shape == (n,)

for scalar-valued Markov parameters or:

data.shape == (n, p, m)

for matrix-valued Markov parameters of dimension \(p\times m\), where \(m\) is the number of inputs and \(p\) is the number of outputs of the system.

sampling_time

A number that denotes the sampling time of the system (in seconds).

force_stability

Whether the Markov parameters are zero-padded to double the length in order to enforce Kung’s stability assumption. See [Kun78]. Defaults to True.

feedthrough

(Optional) Operator or NumPy array of shape (p, m). The zeroth Markov parameter that defines the feedthrough of the realization. Defaults to None.

Methods

input_projector	Construct the right/input projector \(W_2\).
output_projector	Construct the left/output projector \(W_1\).

cache_region = 'memory'

input_projector(num_right)

Construct the right/input projector \(W_2\).

output_projector(num_left)

Construct the left/output projector \(W_1\).

class pymor.reductors.era.RandomizedERAReductor(data, sampling_time, force_stability=True, feedthrough=None, allow_transpose=True, power_iterations=2, rrf_args=None, num_left=None, num_right=None)

Bases: ERAReductorBase

Randomized Eigensystem Realization Algorithm reductor.

Constructs a (reduced) realization from a sequence of Markov parameters \(h_i\), for \(i\in\{1,\,\dots,\,2s-1\}\), \(s\in\mathbb{N}\), by a (reduced) randomized orthogonal factorization of the Hankel matrix of Markov parameters.

See ERAReductor for details on the method and the error bounds. The main difference to ERAReductor is that the conventional SVD of the Hankel matrix is replaced by a randomized SVD for greater efficiency. This implies that any tangential projections for systems with many in- or outputs need to be performed beforehand, i.e. by setting num_left and num_right in the constructor.

The basic randomized ERA algorithm is based on [MSKC21]. Several computational improvements, i.e., a memory-efficient Cholesky QR algorithm [BPS26] and a fast leave-one-out error estimator [ET24], are implemented on top of that. Further, a heuristic error estimator is implemented which enables an adaptive refinement of the ROM that reuses previous computations. The heuristic error estimator, as well as the adaptive numerical algorithm, are described in detail in [PS26].

data

NumPy array that contains the first \(n\) Markov parameters of an LTI system. Has to be one- or three-dimensional with either:

data.shape == (n,)

for scalar-valued Markov parameters or:

data.shape == (n, p, m)

for matrix-valued Markov parameters of dimension \(p\times m\), where \(m\) is the number of inputs and \(p\) is the number of outputs of the system.

sampling_time

A number that denotes the sampling time of the system (in seconds).

force_stability

Whether the Markov parameters are zero-padded to double the length in order to enforce Kung’s stability assumption. See [Kun78]. Defaults to True.

feedthrough

(Optional) Operator or NumPy array of shape (p, m). The zeroth Markov parameter that defines the feedthrough of the realization. Defaults to None.

allow_transpose

Whether to allow the computation of the transposed problem, i.e., the randomized SVD of \(\mathcal{H}^\top\). This can be computationally beneficial if the number of inputs is larger than the number of outputs, as it reduces the second dimension of the Hankel matrix. The transposition only happens internally and is resolved after randomized SVD such that the in- and output dimensionality of the ROM is not affected. Defaults to True.

power_iterations

Number of power iterations used by the randomized SVD. Defaults to 2.

rrf_args

Additional options passed to the underlying randmomized range finder. In particular, qr_method defaults to 'shifted_chol_qr' but can be overridden. The error estimator is always 'loo', as required by the adaptive ERA algorithm, see [PS26].

num_left

Number of left (output) directions for tangential projection.

num_right

Number of right (input) directions for tangential projection.

Methods

reduce	Construct a reduced realization with randomized methods.
relative_error_estimate	Estimate the relative \(\mathcal{H}_2\) error based on the current range space.

reduce(r=None, tol=None)

Construct a reduced realization with randomized methods.

Parameters:

• r – Order of the reduced model if tol is None, maximum order if tol is specified.

• tol – Tolerance for the LOO range-approximation error estimator, scaled by the weighted \(\mathcal{H}_2\) norm. For details, see Section 5.1 in [PS26].

Returns:

rom – Reduced-order LTIModel.

relative_error_estimate()

Estimate the relative \(\mathcal{H}_2\) error based on the current range space.

The estimate is the LOO error estimate for the randomized Hankel approximation, normalized by the weighted \(\mathcal{H}_2\) norm. It is not a rigorous error bound. The latest approximation of the range space will be used for estimation. This means that the error corresponds to a ROM of order r=len(self.randomized_svd.range_finder.Q[-1]). See Section 3.5 in [PS26] for details.

Note

To estimate the error of a ROM whose order is smaller than the latest range-space basis, add the SVD truncation error to this estimate.

Returns:

error – Estimated relative \(\mathcal{H}_2\) error.

Raises:

RuntimeError – If reduce has not been called yet.