SpectralFactorMOR

Represents a system of the form

\[E = \begin{bmatrix} E_{11} & 0 & 0 & 0\\ 0 & E_{22} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad A = \begin{bmatrix} 0 & 0 & 0 & I\\ 0 & A_{22} & 0 & 0\\ 0 & 0 & I & 0\\ -I&0&0&0 \end{bmatrix}\]

with $E_{11}$ and $E_{22}$ nonsingular. The system matrices $B$ and $C$ are partitioned accordingly as

\[B = \begin{bmatrix} B_1\\ B_2\\ B_3\\ B_4 \end{bmatrix}, \quad C = \begin{bmatrix} C_1 & C_2 & C_3 & C_4 \end{bmatrix}.\]

The struct holds options for the IRKA algorithm.

• The conv_tol option sets the convergence tolerance.
• The max_iterations option sets the maximum number of iterations.
• The cycle_detection_length controls how many past iterations are considered for cycle detection. The cycle_detection_tol sets the cycle detection tolerance. If cycle_detection_tol is zero, no cycle detection is performed.
• The s_init_start and s_init_stop control the initial interpolation points. The initial interpolation points are distributed logarithmically between 10^s_init_start and 10^s_init_stop.
• If randomize_s_init is true, the initial interpolation points are disturbed randomly by randn. The variance is set by randomize_s_var.

Represents a system of the form

\[E = \begin{bmatrix} E_{11} & 0\\ 0 & 0 \end{bmatrix}, \quad A = \begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix}\]

with $E_{11}$ and $A_{22}$ nonsingular. The system matrices $B$ and $C$ are partitioned accordingly as

\[B = \begin{bmatrix} B_1\\ B_2 \end{bmatrix}, \quad C = \begin{bmatrix} C_1 & C_2 \end{bmatrix}.\]

Represents a system of the form

\[E = \begin{bmatrix} E_{11} & 0 & 0 & 0\\ 0 & E_{22} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix},\quad A = \begin{bmatrix} A_{11} & A_{12} & A_{13} & A_{14}\\ A_{21} & A_{22} & A_{23} & 0\\ A_{31} & A_{32} & A_{33} & 0\\ A_{41} & 0 & 0 & 0 \end{bmatrix}\]

with $E_{11}, E_{22}, A_{14}=-A_{41}^T,A_{33}$ nonsingular.

Represents a system without holding any structure information (in constract to SemiExplicitIndex1DAE, StaircaseDAE and AlmostKroneckerDAE).

arec_lr_nwt(
    F, E, Q, R, P_r=I, P_l=I; 
    conv_tol=1e-12, max_iterations=20, 
    lyapc_lradi_tol=1e-12, lyapc_lradi_max_iterations=150
)

Solves the Riccati equation

\[FXE^T + EXF^T + EXQ^TQXE^T + P_l RR^T P_l^T = 0,\quad X = P_r X P_r^T.\]

using the low-rank Newton method. Returns an approximate solution $Z$ such that $X ≈ Z Z^T$.

See [1, Alg. 9].

checkprojlure(sys, X)

Checks if X satisfies the projected positive real Lur'e equation, i. e., if there exist matrices $L$ and $M$ such that

\[\begin{align*} \begin{bmatrix} -A^TXE-E^TXA & P_r^TC^T - E^TXB\\ CP_r - B^TXE & M_0+M_0^T \end{bmatrix} &= \begin{bmatrix} L^T\\ M^T \end{bmatrix} \begin{bmatrix} L & M \end{bmatrix},\\ X=X^T &\geq 0,\ X=P_l^TXP_l. \end{align*}\]

compress_lr(Z, r)

Compresses $Z$ via SVD such that it holds $Z Z^T ≈ Z_\mathrm{new} Z_\mathrm{new}^T$, where $Z_\mathrm{new}$ is the returned matrix.

compress_lr(Z; tol=1e-02 * sqrt(size(Z, 1) * eps(Float64)))

Compresses $Z$ via SVD such that it holds $Z Z^T ≈ Z_\mathrm{new} Z_\mathrm{new}^T$, where $Z_\mathrm{new}$ is the returned matrix.

h2normsp(sys::AbstractDescriptorStateSpace)

Computes the H2-norm of the strictly proper system sys.

Uses lyapc_lradi internally.

initial_shifts(A, E, B, max_iterations=10, iter=1)

Returns initial shifts for the LR-ADI method.

See [2, pp. 92-95] and https://github.com/pymor/pymor/blob/main/src/pymor/algorithms/lradi.py.

irka(
    sys::AbstractDescriptorStateSpace, r, irka_options::IRKAOptions;
    P_r=I, P_l=I, Dr=nothing
)

Computes a reduced-order model for sys using the Iterative rational Krylov Algorithm (IRKA) [3].

The parameter r corresponds to the dimension of the reduced-order model. If Dr=nothing, the feed-through matrix $D_r$ of the reduced-order model is set to the original feed-through matrix $D$ of the full-order model. For descriptor systems, the spectral projectors used for interpolation may be provided via the parameters P_r and P_l [4, Sec. 3].

irka(sys::SemiExplicitIndex1DAE, r, irka_options::IRKAOptions)

Computes a reduced-order model for the semi-explicit index-1 system sys using the Iterative rational Krylov Algorithm (IRKA) [3, 4].

The parameter r corresponds to the dimension of the reduced-order model.

irka(sys::StaircaseDAE, r, irka_options::IRKAOptions)

Computes a reduced-order model for the staircase system sys using the Iterative rational Krylov Algorithm (IRKA) [3, 4].

The parameter r corresponds to the dimension of the reduced-order model.

is_valid(sys::StaircaseDAE)

Asserts that sys is in valid staircase form.

isastable(
    sys::Union{SemiExplicitIndex1DAE, AlmostKroneckerDAE, StaircaseDAE}
)

Checks if sys is asymptotically stable, i.e., if all finite eigenvalues lie in the open left-half plane.

ispr(sys::AbstractDescriptorStateSpace) -> (ispr, mineigval, w)

Checks if the system is positive real by sampling the Popov function on the imaginary axis and verifying that the Popov function is positive semi-definite for all sample points.

lrcf(X, trunc_tol)

Computes an approximate low-rank Cholesky-like factorization of a symmetric positive semi-definite matrix $X$ s.t. $X = Z^T Z$ (up to a prescribed tolerance trunc_tol).

lyapc_lradi(A, E, B, P_l=I; tol=1e-12, max_iterations=150)

Computes a low-rank solution for the projected continuous-time Lyapunov equation

\[EXA^T + AXE^T = - P_l BB^T P_l^T,\quad X = P_r X P_r^T\]

using the LR-ADI method [1, Alg. 7]. Returns an approximate solution $Z$ such that $X ≈ Z Z^T$.

new_shifts(A, E, V, Z, prev_shifts, subspace_columns=6)

Returns new intermediate shifts for the LR-ADI method.

See [2, pp. 92-95] and https://github.com/pymor/pymor/blob/main/src/pymor/algorithms/lradi.py.

pr_c_gramian(sys::SemiExplicitIndex1DAE)

Computes the positive real controllability Gramian, i.e., the unique stabilizing solution of

\[AXE^T + EXA^T + (P_l B - EXC^T) (M_0+M_0')^{-1} (P_l B - EXC^T)^T = 0,\quad X = P_r X P_r^T,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Note: Uses the dense solver MatrixEquations.garec.

pr_c_gramian(sys::StaircaseDAE)

Computes the positive real observability Gramian, i.e., the unique stabilizing solution of

\[AXE^T + EXA^T + (P_l B - EXC^T) (M_0+M_0')^{-1} (P_l B - EXC^T)^T = 0,\quad X = P_r X P_r^T,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Note: Uses the dense solver MatrixEquations.garec.

pr_c_gramian_lr(sys::SemiExplicitIndex1DAE)

Computes a low-rank factor of positive real controllability Gramian, i.e., the unique stabilizing solution of

\[AXE^T + EXA^T + (P_l B - EXC^T) (M_0+M_0')^{-1} (P_l B - EXC^T)^T = 0,\quad X = P_r X P_r^T,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Returns an approximate solution $Z$ such that $X ≈ Z Z^T$.

pr_o_gramian(sys::SemiExplicitIndex1DAE)

Computes positive real observability Gramian, i.e., the unique stabilizing solution of

\[A^TXE + E^TXA + (P_r^T C^T-E^TXB)(M_0+M_0')^{-1}(C P_r-B^TXE) = 0,\quad X = P_l^T Y P_l,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Note: Uses the dense solver MatrixEquations.garec.

pr_o_gramian(sys::StaircaseDAE)

Computes positive real observability Gramian, i.e., the unique stabilizing solution of

\[A^TXE + E^TXA + (P_r^T C^T-E^TXB)(M_0+M_0')^{-1}(C P_r-B^TXE) = 0,\quad X = P_l^T Y P_l,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Note: Uses the dense solver MatrixEquations.garec.

pr_o_gramian_lr(sys::SemiExplicitIndex1DAE)

Computes a low-rank factor of the positive real observability Gramian, i.e., the unique stabilizing solution of

\[A^TXE + E^TXA + (P_r^T C^T-E^TXB)(M_0+M_0')^{-1}(C P_r-B^TXE) = 0,\quad X = P_l^T Y P_l,\]

where $P_r$, $P_l$ and $M_0$ are defined by the given system sys.

Returns an approximate solution $Z$ such that $X ≈ Z Z^T$.

prbaltrunc(sys::SemiExplicitIndex1DAE, r, Z_prc, Z_pro)

Computes a reduced-order model for the semi-explicit index-1 system sys using positive real balanced truncation [1, Alg. 2].

The parameter r corresponds to the dimension of the reduced-order model.

The parameter Z_prc must be a (low-rank) Cholesky factor such that $Z_\mathrm{prc}^T Z_\mathrm{prc}$ is the positive real controllability Gramian.

The parameter Z_pro must be a (low-rank) Cholesky factor such that $Z_\mathrm{pro}^T Z_\mathrm{pro}$ is the positive real observability Gramian.

prbaltrunc(sys::StaircaseDAE, r, Z_prc, Z_pro)

Computes a reduced-order model for the staircase system sys using positive real balanced truncation [1, Alg. 2].

The parameter r corresponds to the dimension of the finite part of the reduced-order model.

The parameter Z_prc must be a (low-rank) Cholesky factor such that $Z_\mathrm{prc}^T Z_\mathrm{prc}$ is the positive real controllability Gramian.

The parameter Z_pro must be a (low-rank) Cholesky factor such that $Z_\mathrm{pro}^T Z_\mathrm{pro}$ is the positive real observability Gramian.

retry_irka(run, max_tries, syssp)

Retries IRKA max_tries times and picks the result with the lowest H2-error. The strictly proper subsystem of the full-order model syssp must be supplied in order to compute the H2-error. The parameter run must be a callable which executes IRKA.

sfmor(
    sys::SemiExplicitIndex1DAE, r, irka_options::IRKAOptions;
    X=nothing, compute_factors_tol=1e-12
)

Executes the spectral factor MOR method for a semi-explicit index-1 system.

A solution to the positive real projected Lur'e equation can be supplied via the parameter X. If no X is supplied, the positive real observability Gramian is computed using pr_o_gramian.

sfmor(
    sys::StaircaseDAE, r, irka_options::IRKAOptions;
    X = nothing, compute_factors_tol=1e-12
)

Executes the spectral factor MOR method for a system in staircase form.

A solution to the positive real projected Lur'e equation can be supplied via the parameter X. If no X is supplied, the positive real observability Gramian is computed using pr_o_gramian.

splitsys(sys::SemiExplicitIndex1DAE) -> (sys1, sys2)

Returns the infinite and strictly proper subsystems of sys.

• sys1 is the infinite subsystem.
• sys2 is the strictly proper subsystem.

splitsys(sys::AlmostKroneckerDAE) -> (sys1, sys2)

Returns the infinite and strictly proper subsystems of sys.

• sys1 is the infinite subsystem.
• sys2 is the strictly proper subsystem.

splitsys(sys::StaircaseDAE) -> (sys1, sys2)

Returns the infinite and strictly proper subsystems of sys.

• sys1 is the infinite subsystem.
• sys2 is the strictly proper subsystem.

toindex0(sys::SemiExplicitIndex1DAE)

Transforms sys to a system with index 0. Returns a system of type SemiExplicitIndex1DAE.

tokronecker(sys::StaircaseDAE)

Transforms sys into almost Kronecker form [5, Alg. 5]. Returns a system of type AlmostKroneckerDAE.

The result is cached using the Memoize.jl package.

truncation(d, L, trunc_tol) -> (dr, Lr)

Computes a rank revealing factorization for a given LDL-decomposition of $S = L * \mathrm{diag}(d) * L^T$ (up to a prescribed tolerance trunc_tol) such that $L_r * diag(d_r) * L_r^T \approx S$.