r"""Eigenvalue analysis of the azimuth-averaged state matrix.
Given the non-rotating (MBC-averaged) state matrix :math:`A` in the canonical
``[q2, q2dot, q1]`` layout, an eigenvalue analysis yields the full-system modes:
their natural frequencies, damping ratios, and (complex) mode shapes.
Each under-damped mode appears as a complex-conjugate eigenvalue pair; only the
representative with positive imaginary part is retained. For an eigenvalue
:math:`\lambda`, the modal quantities are
.. math::
\omega_n = |\lambda|, \qquad
\zeta = \frac{-\operatorname{Re}\lambda}{|\lambda|}, \qquad
\omega_d = \operatorname{Im}\lambda,
so the natural and damped frequencies in hertz are
:math:`\omega_n / 2\pi` and :math:`\omega_d / 2\pi`. The velocity (``q2dot``)
rows are dropped from the returned mode shapes, leaving the second-order
displacements and first-order states; the full eigenvectors are also retained.
Notes
-----
OpenFAST exposes no full-system mass or stiffness matrix, so modal frequencies
are obtained directly from :math:`A` and vary with operating point.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import TYPE_CHECKING
import numpy as np
from scipy.linalg import eig as _scipy_eig
if TYPE_CHECKING:
import numpy.typing as npt
from vane.io.mbc_transform import MBCResult
__all__ = ["ModalSolution", "compute_modes", "modes_from_mbc"]
_TWO_PI = 2.0 * np.pi
[docs]
@dataclass
class ModalSolution:
r"""Modal properties extracted from a state matrix.
All per-mode arrays share the same ordering (ascending natural frequency when
the solver sorts). Only under-damped modes (positive imaginary eigenvalue)
are represented.
Parameters
----------
eigenvalues : npt.NDArray[np.complex128]
Retained eigenvalues :math:`\\lambda`, shape ``(n_modes,)``.
natural_frequencies_hz : npt.NDArray[np.float64]
Natural frequencies :math:`|\\lambda| / 2\\pi`, shape ``(n_modes,)``.
damping_ratios : npt.NDArray[np.float64]
Damping ratios :math:`-\\operatorname{Re}\\lambda / |\\lambda|`.
damped_frequencies_hz : npt.NDArray[np.float64]
Damped frequencies :math:`\\operatorname{Im}\\lambda / 2\\pi`.
mode_shapes : npt.NDArray[np.complex128]
Reduced complex mode shapes (second-order displacements then first-order
states), shape ``(ndof2 + ndof1, n_modes)``, phase-normalized so each
mode's dominant component is real and positive.
full_eigenvectors : npt.NDArray[np.complex128]
Full eigenvectors over all states, shape ``(n_states, n_modes)``, sharing
the same phase normalization.
n_rigid_body_modes : int
Number of degrees of freedom without an under-damped representative
(rigid-body or over-damped modes).
dof_descriptions : list[str]
Descriptions of the degrees of freedom that index the rows of
``mode_shapes`` (the second-order displacements then first-order states).
Empty when descriptions were not supplied to the solver.
dof_mbc_coordinates : list[str]
Multi-blade coordinate (``"collective"``/``"cosine"``/``"sine"``/``""``) of
each ``mode_shapes`` row, aligned with ``dof_descriptions``. Empty when not
supplied.
condition_numbers : npt.NDArray[np.float64]
Per-mode eigenvalue condition number :math:`\\kappa = 1 / |y^H x|` (the
reciprocal cosine of the angle between the left and right eigenvectors),
shape ``(n_modes,)``. ``1`` is perfectly conditioned; large values flag
eigenvalues that are numerically sensitive (highly non-normal system).
is_degenerate : npt.NDArray[np.bool_]
Per-mode flag marking modes whose eigenvalue is (near-)repeated by another
retained mode, so the individual mode shape is defined only up to a
rotation within the shared invariant subspace.
n_unstable : int
Number of eigenvalues with positive real part (unstable). Non-zero is a
red flag for the operating point.
n_overdamped : int
Number of non-oscillatory stable real eigenvalues (over-damped modes),
which are not returned as oscillatory modes.
"""
eigenvalues: npt.NDArray[np.complex128]
natural_frequencies_hz: npt.NDArray[np.float64]
damping_ratios: npt.NDArray[np.float64]
damped_frequencies_hz: npt.NDArray[np.float64]
mode_shapes: npt.NDArray[np.complex128]
full_eigenvectors: npt.NDArray[np.complex128]
n_rigid_body_modes: int
dof_descriptions: list[str]
condition_numbers: npt.NDArray[np.float64] = field(
default_factory=lambda: np.empty(0, dtype=np.float64)
)
is_degenerate: npt.NDArray[np.bool_] = field(
default_factory=lambda: np.empty(0, dtype=np.bool_)
)
n_unstable: int = 0
n_overdamped: int = 0
# Appended after the Stage-1 optional fields to keep positional construction
# of the previously-existing fields stable.
dof_mbc_coordinates: list[str] = field(default_factory=list)
@property
def n_modes(self) -> int:
"""Return the number of retained modes."""
return int(self.eigenvalues.shape[0])
@property
def natural_frequencies_rad_s(self) -> npt.NDArray[np.float64]:
"""Return natural frequencies in rad/s."""
return self.natural_frequencies_hz * _TWO_PI
@property
def max_condition_number(self) -> float:
"""Return the worst (largest) eigenvalue condition number, or 1 if none."""
if self.condition_numbers.size == 0:
return 1.0
return float(self.condition_numbers.max())
[docs]
def compute_modes(
a: npt.NDArray[np.float64],
ndof2: int,
ndof1: int,
*,
sort_by_frequency: bool = True,
descriptions: list[str] | None = None,
mbc_coordinates: list[str] | None = None,
) -> ModalSolution:
"""Compute modal properties from a state matrix.
Parameters
----------
a : npt.NDArray[np.float64]
Square state matrix in ``[q2, q2dot, q1]`` order, shape
``(2 * ndof2 + ndof1, 2 * ndof2 + ndof1)``.
ndof2 : int
Number of second-order degrees of freedom (the ``q2`` block size).
ndof1 : int
Number of first-order states (the ``q1`` block size).
sort_by_frequency : bool, optional
Sort the returned modes by ascending natural frequency (the default).
descriptions : list[str] or None, optional
Descriptions of all ``2 * ndof2 + ndof1`` states, in the matrix order.
When given, the subset that indexes the mode-shape rows is stored on the
result as ``dof_descriptions``.
mbc_coordinates : list[str] or None, optional
Multi-blade coordinate of all states, in the matrix order; the mode-shape
subset is stored on the result as ``dof_mbc_coordinates``.
Returns
-------
ModalSolution
The modal frequencies, damping ratios, and mode shapes.
Raises
------
ValueError
If ``a`` is not square or its size does not match ``2 * ndof2 + ndof1``.
"""
if a.ndim != 2 or a.shape[0] != a.shape[1]:
msg = f"State matrix must be square, got shape {a.shape}"
raise ValueError(msg)
n_states = a.shape[0]
if 2 * ndof2 + ndof1 != n_states:
msg = (
f"State count mismatch: 2*ndof2 + ndof1 = {2 * ndof2 + ndof1} "
f"but A is {n_states}x{n_states}"
)
raise ValueError(msg)
if not np.all(np.isfinite(a)):
msg = "State matrix contains non-finite values (NaN or Inf)"
raise ValueError(msg)
eigen = _scipy_eig(a, left=True, right=True)
eigenvalues = np.asarray(eigen[0], dtype=np.complex128)
left_vectors = np.asarray(eigen[1], dtype=np.complex128)
right_vectors = np.asarray(eigen[2], dtype=np.complex128)
n_unstable, n_overdamped = _classify_spectrum(eigenvalues)
conditioning = _eigenvalue_conditioning(left_vectors, right_vectors)
keep = np.flatnonzero(eigenvalues.imag > 0.0)
lam = eigenvalues[keep].astype(np.complex128)
condition = conditioning[keep]
degenerate = _degeneracy_flags(lam)
omega_n = np.abs(lam)
damping = -lam.real / omega_n
omega_d = lam.imag
displacement_rows = np.r_[0:ndof2, 2 * ndof2 : n_states]
reduced = right_vectors[np.ix_(displacement_rows, keep)].astype(np.complex128)
full = right_vectors[:, keep].astype(np.complex128)
_phase_normalize(reduced, full)
if sort_by_frequency:
order = np.argsort(omega_n, kind="stable")
lam, omega_n, damping, omega_d, condition, degenerate = (
lam[order],
omega_n[order],
damping[order],
omega_d[order],
condition[order],
degenerate[order],
)
reduced, full = reduced[:, order], full[:, order]
reduced_rows = displacement_rows.tolist()
dof_descriptions = (
[descriptions[i] for i in reduced_rows] if descriptions is not None else []
)
dof_mbc_coordinates = (
[mbc_coordinates[i] for i in reduced_rows] if mbc_coordinates else []
)
return ModalSolution(
eigenvalues=lam,
natural_frequencies_hz=omega_n / _TWO_PI,
damping_ratios=damping,
damped_frequencies_hz=omega_d / _TWO_PI,
mode_shapes=reduced,
full_eigenvectors=full,
n_rigid_body_modes=(ndof2 + ndof1) - int(keep.shape[0]),
dof_descriptions=dof_descriptions,
dof_mbc_coordinates=dof_mbc_coordinates,
condition_numbers=condition,
is_degenerate=degenerate,
n_unstable=n_unstable,
n_overdamped=n_overdamped,
)
[docs]
def modes_from_mbc(result: MBCResult) -> ModalSolution:
"""Compute modal properties from an MBC transform result.
Parameters
----------
result : MBCResult
Output of :func:`vane.io.mbc_transform.mbc3_transform`; its ``avg_a`` must
be present.
Returns
-------
ModalSolution
The modal solution for the averaged state matrix.
Raises
------
ValueError
If ``result.avg_a`` is ``None``.
"""
if result.avg_a is None:
msg = "MBCResult has no averaged A matrix to analyse"
raise ValueError(msg)
return compute_modes(
result.avg_a,
result.ndof2,
result.ndof1,
descriptions=result.state_descriptions,
mbc_coordinates=result.mbc_coordinates or None,
)
def _phase_normalize(
reduced: npt.NDArray[np.complex128], full: npt.NDArray[np.complex128]
) -> None:
"""Rotate each mode's phase in place so its dominant component is real positive.
For every column, the component of largest magnitude in ``reduced`` is made
real and positive by dividing the column (in both ``reduced`` and ``full``) by
that component's unit phase. This gives a deterministic, sign-stable mode
shape across operating points.
Parameters
----------
reduced : npt.NDArray[np.complex128]
Reduced mode-shape matrix, modified in place.
full : npt.NDArray[np.complex128]
Full eigenvector matrix, modified in place with the same rotation.
"""
for col in range(reduced.shape[1]):
dominant = reduced[np.argmax(np.abs(reduced[:, col])), col]
magnitude = np.abs(dominant)
if magnitude == 0.0:
continue
unit_phase = dominant / magnitude
reduced[:, col] /= unit_phase
full[:, col] /= unit_phase
def _classify_spectrum(eigenvalues: npt.NDArray[np.complex128]) -> tuple[int, int]:
"""Count unstable and over-damped eigenvalues in the full spectrum.
Parameters
----------
eigenvalues : npt.NDArray[np.complex128]
All eigenvalues of the state matrix.
Returns
-------
tuple[int, int]
``(n_unstable, n_overdamped)`` — eigenvalues with positive real part, and
non-zero stable real (non-oscillatory) eigenvalues, using a tolerance
scaled to the spectral radius.
"""
if eigenvalues.size == 0:
return 0, 0
# Absolute numerical-zero floor (~1000x the eigenvalue round-off, eps*scale),
# NOT a coarse fraction of the spectral radius: scaling the threshold by the
# largest eigenvalue would hide slow instabilities in a stiff spectrum.
scale = max(float(np.abs(eigenvalues).max()), 1.0)
tolerance = 1e3 * float(np.finfo(np.float64).eps) * scale
n_unstable = int(np.count_nonzero(eigenvalues.real > tolerance))
is_real = np.abs(eigenvalues.imag) <= tolerance
is_nonzero = np.abs(eigenvalues) > tolerance
n_overdamped = int(
np.count_nonzero(is_real & (eigenvalues.real < -tolerance) & is_nonzero)
)
return n_unstable, n_overdamped
def _eigenvalue_conditioning(
left: npt.NDArray[np.complex128], right: npt.NDArray[np.complex128]
) -> npt.NDArray[np.float64]:
r"""Return each eigenvalue's condition number from its left/right eigenvectors.
The condition number is :math:`\kappa = 1 / |y^H x|` with unit-norm left
eigenvector ``y`` and right eigenvector ``x``; by Cauchy-Schwarz ``kappa >= 1``,
equal to 1 for a normal matrix and large for an ill-conditioned (highly
non-normal) eigenvalue.
Parameters
----------
left, right : npt.NDArray[np.complex128]
Unit-norm left and right eigenvectors as columns, in matching order.
Returns
-------
npt.NDArray[np.float64]
The condition number per eigenvalue (``inf`` for a defective eigenvalue).
"""
sensitivity = np.abs(np.sum(np.conj(left) * right, axis=0))
with np.errstate(divide="ignore"):
condition = np.where(sensitivity > 0.0, 1.0 / sensitivity, np.inf)
return np.asarray(condition, dtype=np.float64)
def _degeneracy_flags(
eigenvalues: npt.NDArray[np.complex128], *, rtol: float = 1e-6
) -> npt.NDArray[np.bool_]:
"""Flag eigenvalues that are (near-)repeated by another in the set.
Parameters
----------
eigenvalues : npt.NDArray[np.complex128]
The retained eigenvalues.
rtol : float, optional
Relative tolerance for treating two eigenvalues as coincident.
Returns
-------
npt.NDArray[np.bool_]
``True`` for each eigenvalue within ``rtol`` of another (its mode shape is
then defined only up to a rotation within the shared invariant subspace).
"""
n = eigenvalues.shape[0]
flags = np.zeros(n, dtype=np.bool_)
for i in range(n):
for j in range(i + 1, n):
scale = max(abs(eigenvalues[i]), abs(eigenvalues[j]), 1e-30)
if abs(eigenvalues[i] - eigenvalues[j]) <= rtol * scale:
flags[i] = True
flags[j] = True
return flags