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.StateSpace) – State space to reduce

Returns

Reduced state space system

Return type

(libss.StateSpace)

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

save_reduced_order_bases(file_name)[source]

Save reduced order bases to an h5 file

Parameters

file_name (str) – path to h5 file

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.

two_sided_arnoldi(frequency, r)[source]

Two-sided projection with a single interpolation point following the Arnoldi procedure. Very similar to the one-sided method available, but it adds the projection \(\mathbf{W}\) built using the Krylov space for the \(\mathbf{c}\) vector:

\[\mathcal{K}_r((\sigma\mathbf{I}_n - \mathbf{A})^{-T}, (\sigma\mathbf{I}_n - \mathbf{A})^{-T}\mathbf{c}^T)\subseteq\mathcal{W}=\text{range}(\mathbf{W})\]

The oblique projection \(\mathbf{VW}^T\) matches twice as many moments as the single sided projection.

The resulting system takes the form:

\[\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}=W^TAV\in\mathbb{R}^{k\times k},\,\\ \hat{B}=W^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 on each side. The resulting ROM will be of order \(2r\).

Returns

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

Return type

tuple