Krylov¶

class sharpy.rom.krylov.Krylov[source]¶

Model Order Reduction Methods for Single Input Single Output (SISO) and MIMO Linear Time-Invariant (LTI) Systems using moment matching (Krylov Methods).

Examples

General calling sequences for different systems

SISO single point interpolation:

>>> algorithm = 'one_sided_arnoldi'
>>> interpolation_point = np.array([0.0])
>>> krylov_r = 4
>>>
>>> rom = Krylov()
>>> rom.initialise(sharpy_data, FullOrderModelSS)
>>> rom.run(algorithm, krylov_r, interpolation_point)

2 by 2 MIMO with tangential, multipoint interpolation:

>>> algorithm = 'dual_rational_arnoldi'
>>> interpolation_point = np.array([0.0, 1.0j])
>>> krylov_r = 4
>>> right_vector = np.block([[1, 0], [0, 1]])
>>> left_vector = right_vector
>>>
>>> rom = Krylov()
>>> rom.initialise(sharpy_data, FullOrderModelSS)
>>> rom.run(algorithm, krylov_r, interpolation_point, right_vector, left_vector)

2 by 2 MIMO multipoint interpolation:

>>> algorithm = 'mimo_rational_arnoldi'
>>> interpolation_point = np.array([0.0])
>>> krylov_r = 4
>>>
>>> rom = Krylov()
>>> rom.initialise(sharpy_data, FullOrderModelSS)
>>> rom.run(algorithm, krylov_r, interpolation_point)

The settings that this solver accepts are given by a dictionary, with the following key-value pairs:

Name	Type	Description	Default	Options
`print_info`	`bool`	Write ROM information to screen and log	`True`
`frequency`	`list(complex)`	Interpolation points in the continuous time complex plane [rad/s]	`[0]`
`algorithm`	`str`	Krylov reduction method algorithm
`r`	`int`	Moments to match at the interpolation points	`1`
`single_side`	`str`	Construct the rom using a single side. Leave blank (or empty string) for both.		`controllability`, `observability`
`tangent_input_file`	`str`	Filepath to .h5 file containing tangent interpolation vectors
`restart_arnoldi`	`bool`	Restart Arnoldi iteration with r-=1 if ROM is unstable	`False`

check_stability(restart_arnoldi=False)[source]¶

Checks the stability of the ROM by computing its eigenvalues.

If the resulting system is unstable, the Arnoldi procedure can be restarted to eliminate the eigenvalues outside the stability boundary.

However, if this is the case, the ROM no longer matches the moments of the original system at the specific frequencies since now the approximation is done with respect to a system of the form:

\[\begin{split}\Sigma = \left(\begin{array}{c|c} \mathbf{A} & \mathbf{\bar{B}} \\ \hline \mathbf{C} & \ \end{array}\right)\end{split}\]

where \(\mathbf{\bar{B}} = (\mu \mathbf{I}_n - \mathbf{A})\mathbf{B}\)

Parameters: restart_arnoldi (bool) – Restart the relevant Arnoldi algorithm with the unstable eigenvalues removed.

dual_rational_arnoldi(frequency, r)[source]¶

Dual Rational Arnoli Interpolation for SISO sytems [1] and MIMO systems through tangential interpolation [2].

Effectively the same as the two_sided_arnoldi and the resulting V matrices for each interpolation point are concatenated

\[\begin{split}\bigcup\limits_{k = 1}^K\mathcal{K}_{b_k}((\sigma_i\mathbf{I}_n - \mathbf{A})^{-1}, (\sigma_i\mathbf{I}_n - \mathbf{A})^{-1}\mathbf{b})\subseteq\mathcal{V}&=\text{range}(\mathbf{V}) \\ \bigcup\limits_{k = 1}^K\mathcal{K}_{c_k}((\sigma_i\mathbf{I}_n - \mathbf{A})^{-T}, (\sigma_i\mathbf{I}_n - \mathbf{A})^{-T}\mathbf{c}^T)\subseteq\mathcal{Z}&=\text{range}(\mathbf{Z})\end{split}\]

For MIMO systems, tangential interpolation is used through the right and left tangential direction vectors \(\mathbf{r}_i\) and \(\mathbf{l}_i\).

\[\begin{split}\bigcup\limits_{k = 1}^K\mathcal{K}_{b_k}((\sigma_i\mathbf{I}_n - \mathbf{A})^{-1}, (\sigma_i\mathbf{I}_n - \mathbf{A})^{-1}\mathbf{Br}_i)\subseteq\mathcal{V}&=\text{range}(\mathbf{V}) \\ \bigcup\limits_{k = 1}^K\mathcal{K}_{c_k}((\sigma_i\mathbf{I}_n - \mathbf{A})^{-T}, (\sigma_i\mathbf{I}_n - \mathbf{A})^{-T}\mathbf{C}^T\mathbf{l}_i)\subseteq\mathcal{Z}&=\text{range}(\mathbf{Z})\end{split}\]

Parameters

frequency (np.ndarray) – Array containing the interpolation points \(\sigma = \{\sigma_1, \dots, \sigma_K\}\in\mathbb{C}\)
r (int) – Krylov space order \(b_k\) and \(c_k\). At the moment, different orders for the controllability and observability constructions are not supported.
right_tangent (np.ndarray) – Matrix containing the right tangential direction interpolation vector for each interpolation point in column form, i.e. \(\mathbf{r}\in\mathbb{R}^{m \times K}\).
left_tangent (np.ndarray) – Matrix containing the left tangential direction interpolation vector for each interpolation point in column form, i.e. \(\mathbf{l}\in\mathbb{R}^{p \times K}\).

Returns

The reduced order model matrices: \(\mathbf{A}_r\), \(\mathbf{B}_r\) and \(\mathbf{C}_r\).

Return type

tuple

References

[1] Grimme [2] Gallivan

mimo_rational_arnoldi(frequency, r)[source]¶

Construct full rank orthonormal projection basis \(\mathbf{V}\) and \(\mathbf{W}\).

The main issue that one normally encounters with MIMO systems is that the minimality assumption of the system does not guarantee the resulting Krylov space to be full rank, unlike in the SISO case. Therefore, the construction is performed vector by vector, where linearly dependent vectors are eliminated or deflated from the Krylov subspace.

If the number of inputs differs the number of outputs, both Krylov spaces will be built such that both are the same size, therefore one Krylov space may be of higher order than the other one.

Following the method for vector-wise construction in Gugercin [1].

Parameters

frequency (np.ndarray) – Array containing interpolation frequencies
r (int) – Krylov space order

Returns

Tuple of reduced system matrices A, B and C.

Return type

tuple

References

[1] Gugercin, S. Projection Methods for Model Reduction of Large-Scale Dynamical: Systems PhD Thesis. Rice University 2003.

one_sided_arnoldi(frequency, r)[source]¶

One-sided Arnoldi method expansion about a single interpolation point, \(\sigma\). The projection matrix \(\mathbf{V}\) is constructed using an order \(r\) Krylov space. The space for a single finite interpolation point known as a Pade approximation is described by:

\[\text{range}(\textbf{V}) = \mathcal{K}_r((\sigma\mathbf{I}_n - \mathbf{A})^{-1}, (\sigma\mathbf{I}_n - \mathbf{A})^{-1}\mathbf{b})\]

In the case of an interpolation about infinity, the problem is known as partial realisation and the Krylov space is

\[\text{range}(\textbf{V}) = \mathcal{K}_r(\mathbf{A}, \mathbf{b})\]

The resulting orthogonal projection leads to the following reduced order system:

\[\begin{split}\hat{\Sigma} : \left(\begin{array}{c|c} \hat{A} & \hat{B} \\ \hline \hat{C} & {D}\end{array}\right) \text{with } \begin{cases}\hat{A}=V^TAV\in\mathbb{R}^{k\times k},\,\\ \hat{B}=V^TB\in\mathbb{R}^{k\times m},\,\\ \hat{C}=CV\in\mathbb{R}^{p\times k},\,\\ \hat{D}=D\in\mathbb{R}^{p\times m}\end{cases}\end{split}\]

Parameters

frequency (complex) – Interpolation point \(\sigma \in \mathbb{C}\)
r (int) – Number of moments to match. Equivalent to Krylov space order and order of the ROM.

Returns

The reduced order model matrices: \(\mathbf{A}_r\), \(\mathbf{B}_r\) and \(\mathbf{C}_r\)

Return type

tuple

real_rational_arnoldi(frequency, r)[source]¶

When employing complex frequencies, the projection matrix can be normalised to be real Following Algorithm 1b in Lee(2006) :param frequency: :param r:

Returns:

restart()[source]¶: Implicitly Restarted Krylov Algorithm

run(ss)[source]¶

Performs Model Order Reduction employing Krylov space projection methods.

Supported methods include:

Algorithm	Interpolation Points	Systems
`one_sided_arnoldi`	1	SISO Systems
`two_sided_arnoldi`	1	SISO Systems
`dual_rational_arnoldi`	K	SISO systems and Tangential interpolation for MIMO systems
`mimo_rational_arnoldi`	K	MIMO systems. Uses vector-wise construction (more robust)
`mimo_block_arnoldi`	K	MIMO systems. Uses block Arnoldi methods (more efficient)

Parameters: ss (sharpy.linear.src.libss.ss) – State space to reduce
Returns: Reduced state space system
Return type: (libss.ss)

save(filename)[source]¶

Saves to an .h5 file of name filename the left and right projectors and the reduced order model

Parameters: filename (str) – path and filename to which to save the data

stable_realisation(*args, **kwargs)[source]¶

Remove unstable poles left after reduction

Using a Schur decomposition of the reduced plant matrix \(\mathbf{A}_m\in\mathbb{C}^{m\times m}\), the method removes the unstable eigenvalues that could have appeared after the moment-matching reduction.

The oblique projection matrices \(\mathbf{T}_L\in\mathbb{C}^{m \times p}\) and \(\mathbf{T}_R\in\mathbb{C}^{m \times p}\) result in a stable realisation

\[\mathbf{A}_s = \mathbf{T}_L^\top\mathbf{AT}_R \in \mathbb{C}^{p\times p}.\]

Parameters

A (np.ndarray) – plant matrix (if not provided self.ssrom.A will be used).

Returns

Left and right projection matrices \(\mathbf{T}_L\in\mathbb{C}^{m \times p}\) and: \(\mathbf{T}_R\in\mathbb{C}^{m \times p}\)

Return type

tuple

References

Jaimoukha, I. M., Kasenally, E. D.. Implicitly Restarted Krylov Subspace Methods for Stable Partial Realizations. SIAM Journal of Matrix Analysis and Applications, 1997.