Inter-Area Oscillation Modes — Matrix-Pencil Damping Estimator

monitor.oscillation_modes recovers the electromechanical modes of a grid ringdown — the response of a wide-area signal (a bus frequency, a tie-line angle, the network order parameter) after a disturbance. The ringdown is a sum of damped sinusoids; each is a mode, and a mode's damping ratio is the quantity that matters for reliability: NERC PRC-028-1 disturbance monitoring and PRC-030-1 unexpected IBR event-mitigation evidence workflows depend on measured disturbance records and post-event analysis, so this estimator reports the frequency and damping evidence that reviewer packages consume.

Method

estimate_oscillation_modes uses the matrix-pencil method of Hua & Sarkar (1990):

1. Build a Hankel matrix from the ringdown samples.
2. Take its singular value decomposition (LAPACK via NumPy) and keep the dominant signal subspace — the model order, chosen from the singular spectrum or set explicitly.
3. Recover each discrete pole z = exp((−α + j·2πf)·Δt) as a generalised eigenvalue of the pencil formed from the right singular vectors.
4. Map each pole to a physical mode: frequency = |∠z|·fs / 2π, ζ = −ln|z| / hypot(ln|z|, ∠z), and recover amplitude and phase from a least-squares Vandermonde fit.

Matrix pencil is a one-shot SVD problem rather than a polynomial root-find, so it is far less noise-sensitive than Prony. Complex poles of a real signal occur in conjugate pairs, which are merged into one positive-frequency mode; a real pole is reported as a pure decay at frequency = 0. Modes are returned ordered by descending amplitude, and any mode whose damping ratio is below the screening threshold (DEFAULT_DAMPING_THRESHOLD = 0.03) is flagged poorly_damped; a growing (unstable) mode has a negative damping ratio.

Evidence boundary

The public estimator validates the ringdown before NumPy conversion: samples must be a one-dimensional finite real sequence with at least four values. Boolean, complex, numeric-text, coercive object, overflowing, and broken array- protocol inputs fail closed with ValueError; legitimate Python and NumPy real numeric object scalars remain supported.

OscillationMode applies the same custody rule to direct construction, not only to records emitted by the estimator. Frequency and amplitude must be finite and non-negative, damping ratio must lie in [-1, 1], phase must lie in [-π, π], and poorly_damped must be a canonical boolean. Accepted real scalar types are normalised to plain JSON-safe float values before to_dict() can publish the record. The flag itself is threshold-dependent, so a directly constructed record can validate its type but cannot replay a caller-specific threshold that is not stored in the record.

Mode-family Screen

Each OscillationMode.to_dict() record now carries mode_family, computed by classify_oscillation_band. The default taxonomy is an engineering review screen, not a regulatory clause threshold:

• aperiodic — near-zero-frequency real-pole decay;
• inter_area — low-frequency area-to-area swings below 1 Hz;
• local — local electromechanical swings below 3 Hz;
• sub_synchronous — higher oscillations below the configured synchronous grid frequency (60 Hz by default);
• super_synchronous — modes at or above the configured synchronous frequency.

The cut-points are explicit function parameters so 50 Hz systems or operator practice can override them without changing the estimator. The PRC evidence record preserves per-mode families and aggregates mode_family_counts, allowing one package to distinguish inter-area and sub-synchronous review signals while keeping the same review-only claim boundary.

Relationship to autotune.freq_id

autotune.freq_id runs multichannel DMD to assign oscillator channels to modal frequencies for tuning. This estimator answers a different question for oscillation safety monitoring: the modal damping of a single ringdown signal. Both are one-shot offline SVD/eigen-solves on a short window (the NumPy floor — not a per-step hot path), so neither carries the multi-language acceleration chain.

Review-only

Like every monitor primitive, the estimator only reads a signal and reports modes; it never changes bindings, layers, or coupling.

oscillation_modes

Inter-area oscillation mode estimation by the matrix-pencil method.

After a grid disturbance, the ringdown of a wide-area signal (a bus frequency, a tie-line angle, the network order parameter) is a sum of damped sinusoids; each sinusoid is an electromechanical mode. A mode's damping ratio is the reliability quantity that matters — North American reliability standards (NERC PRC-028-1 / PRC-030-1 evidence workflows) treat a mode whose damping ratio falls below a few percent as a poorly-damped, flaggable oscillation.

:func:estimate_oscillation_modes fits the ringdown with the matrix-pencil method of Hua & Sarkar (1990): build a Hankel matrix from the samples, take its singular value decomposition, keep the dominant subspace, and recover each mode's discrete pole z = exp((−α + j·2πf)·Δt) as a generalised eigenvalue of the pencil formed from the right singular vectors. The pole gives the modal frequency f = angle(z)·fs / 2π and damping ratio ζ = −ln|z| / hypot(ln|z|, angle(z)); a least-squares Vandermonde fit gives each mode's amplitude and phase. The estimator is diagnostic only — it reads a signal and reports modes; it never changes bindings, layers, or coupling.

Matrix pencil is preferred over Prony for ringdown analysis (it is a one-shot SVD problem rather than a polynomial root-find, far less noise-sensitive) and is distinct in purpose from autotune.freq_id (multichannel DMD that assigns oscillator channels to modal frequencies for tuning): this estimator reports single-signal modal damping for oscillation safety monitoring.

The whole estimate is one offline SVD + eigen-solve on a short ringdown window (LAPACK via NumPy), not a per-step hot path, so it stays on the NumPy floor — the same judgement as autotune.freq_id.

References

• Hua, Y. & Sarkar, T. K. 1990, IEEE Trans. Acoust. Speech Signal Process. 38(5):814–824 — matrix pencil for estimating parameters of exponentially damped/undamped sinusoids in noise.
• Sarkar, T. K. & Pereira, O. 1995, IEEE Antennas Propag. Mag. 37(1):48–55 — using the matrix-pencil method to estimate the parameters of a sum of complex exponentials.
• NERC PRC-028 (oscillation monitoring) — damping-ratio screening of inter-area modes.

Classes

OscillationMode dataclass

OscillationMode(
    frequency_hz: float,
    damping_ratio: float,
    amplitude: float,
    phase_rad: float,
    poorly_damped: bool,
)

A single damped-sinusoid mode recovered from a ringdown.

Attributes

frequency_hz : float Modal oscillation frequency in hertz (≥ 0). damping_ratio : float Dimensionless damping ratio ζ; > 0 is stable, ≤ 0 is growing (unstable), small positive values are poorly damped. amplitude : float Modal amplitude in the units of the input signal (≥ 0). phase_rad : float Modal phase at the first sample, in radians on (−π, π]. poorly_damped : bool Whether damping_ratio is below the screening threshold.

Methods:

__post_init__

__post_init__() -> None

Validate and normalise directly constructed modal evidence.

Source code in src/scpn_phase_orchestrator/monitor/oscillation_modes.py

def __post_init__(self) -> None:
    """Validate and normalise directly constructed modal evidence."""
    frequency = _non_negative_real(self.frequency_hz, "frequency_hz")
    damping = _real_scalar(self.damping_ratio, "damping_ratio")
    if not -1.0 <= damping <= 1.0:
        raise ValueError("damping_ratio must lie in [-1, 1]")
    amplitude = _non_negative_real(self.amplitude, "amplitude")
    phase = _real_scalar(self.phase_rad, "phase_rad")
    if not -np.pi <= phase <= np.pi:
        raise ValueError("phase_rad must lie in [-pi, pi]")
    poorly_damped = _plain_bool(self.poorly_damped, "poorly_damped")
    object.__setattr__(self, "frequency_hz", frequency)
    object.__setattr__(self, "damping_ratio", damping)
    object.__setattr__(self, "amplitude", amplitude)
    object.__setattr__(self, "phase_rad", phase)
    object.__setattr__(self, "poorly_damped", poorly_damped)

to_dict

to_dict() -> dict[str, bool | float | str]

Return a JSON-serialisable mapping of the mode.

Returns

dict[str, bool | float | str] The frequency, damping ratio, amplitude, phase, and the poorly-damped flag plus the engineering mode family.

Source code in src/scpn_phase_orchestrator/monitor/oscillation_modes.py

def to_dict(self) -> dict[str, bool | float | str]:
    """Return a JSON-serialisable mapping of the mode.

    Returns
    -------
    dict[str, bool | float | str]
        The frequency, damping ratio, amplitude, phase, and the
        poorly-damped flag plus the engineering mode family.
    """
    return {
        "frequency_hz": self.frequency_hz,
        "damping_ratio": self.damping_ratio,
        "amplitude": self.amplitude,
        "phase_rad": self.phase_rad,
        "poorly_damped": self.poorly_damped,
        "mode_family": classify_oscillation_band(self.frequency_hz),
    }

Functions:

classify_oscillation_band

classify_oscillation_band(
    frequency_hz: float,
    *,
    inter_area_max_hz: float = 1.0,
    local_max_hz: float = 3.0,
    synchronous_frequency_hz: float = 60.0,
    zero_frequency_tolerance_hz: float = 1e-09,
) -> str

Classify a detected mode into a grid oscillation family.

The defaults are an engineering review taxonomy, not regulatory clause thresholds: near-zero modes are aperiodic decay, modes below 1 Hz are inter-area swings, modes below 3 Hz are local electromechanical swings, and higher modes below the configured synchronous frequency are labelled sub-synchronous for review. Callers may override the cut-points for 50 Hz systems or local operator practice.

Parameters

frequency_hz : float Non-negative modal frequency in hertz. inter_area_max_hz : float Exclusive upper bound for the inter-area family. local_max_hz : float Exclusive upper bound for the local electromechanical family. synchronous_frequency_hz : float Nominal grid synchronous frequency; modes below it and above local_max_hz are labelled sub-synchronous. zero_frequency_tolerance_hz : float Non-negative tolerance below which a mode is treated as aperiodic.

Returns

str One of :data:APERIODIC_MODE, :data:INTER_AREA_MODE, :data:LOCAL_MODE, :data:SUB_SYNCHRONOUS_MODE, or :data:SUPER_SYNCHRONOUS_MODE.

Raises

ValueError If any frequency or boundary is non-finite, negative, or unordered.

Source code in src/scpn_phase_orchestrator/monitor/oscillation_modes.py

def classify_oscillation_band(
    frequency_hz: float,
    *,
    inter_area_max_hz: float = 1.0,
    local_max_hz: float = 3.0,
    synchronous_frequency_hz: float = 60.0,
    zero_frequency_tolerance_hz: float = 1.0e-9,
) -> str:
    """Classify a detected mode into a grid oscillation family.

    The defaults are an engineering review taxonomy, not regulatory clause
    thresholds: near-zero modes are aperiodic decay, modes below 1 Hz are
    inter-area swings, modes below 3 Hz are local electromechanical swings,
    and higher modes below the configured synchronous frequency are labelled
    sub-synchronous for review. Callers may override the cut-points for
    50 Hz systems or local operator practice.

    Parameters
    ----------
    frequency_hz : float
        Non-negative modal frequency in hertz.
    inter_area_max_hz : float
        Exclusive upper bound for the inter-area family.
    local_max_hz : float
        Exclusive upper bound for the local electromechanical family.
    synchronous_frequency_hz : float
        Nominal grid synchronous frequency; modes below it and above
        ``local_max_hz`` are labelled sub-synchronous.
    zero_frequency_tolerance_hz : float
        Non-negative tolerance below which a mode is treated as aperiodic.

    Returns
    -------
    str
        One of :data:`APERIODIC_MODE`, :data:`INTER_AREA_MODE`,
        :data:`LOCAL_MODE`, :data:`SUB_SYNCHRONOUS_MODE`, or
        :data:`SUPER_SYNCHRONOUS_MODE`.

    Raises
    ------
    ValueError
        If any frequency or boundary is non-finite, negative, or unordered.
    """
    frequency = _non_negative_real(frequency_hz, "frequency_hz")
    inter_area_max = _positive_real(inter_area_max_hz, "inter_area_max_hz")
    local_max = _positive_real(local_max_hz, "local_max_hz")
    synchronous_frequency = _positive_real(
        synchronous_frequency_hz, "synchronous_frequency_hz"
    )
    zero_tolerance = _non_negative_real(
        zero_frequency_tolerance_hz, "zero_frequency_tolerance_hz"
    )
    if not zero_tolerance < inter_area_max < local_max < synchronous_frequency:
        raise ValueError(
            "frequency bands must satisfy zero_tolerance < inter_area_max_hz "
            "< local_max_hz < synchronous_frequency_hz"
        )
    if frequency <= zero_tolerance:
        return APERIODIC_MODE
    if frequency < inter_area_max:
        return INTER_AREA_MODE
    if frequency < local_max:
        return LOCAL_MODE
    if frequency < synchronous_frequency:
        return SUB_SYNCHRONOUS_MODE
    return SUPER_SYNCHRONOUS_MODE

estimate_oscillation_modes

estimate_oscillation_modes(
    signal: FloatArray,
    fs: float,
    *,
    model_order: int | None = None,
    pencil_factor: float = 0.4,
    damping_threshold: float = DEFAULT_DAMPING_THRESHOLD,
    energy_floor: float = 0.001,
) -> tuple[OscillationMode, ...]

Estimate damped oscillation modes of a ringdown by the matrix-pencil method.

Parameters

signal : FloatArray Real ringdown samples, uniformly sampled, length ≥ 4. fs : float Sampling frequency in hertz (> 0). model_order : int | None Number of modes (real exponentials plus conjugate pairs) to recover. When None the order is chosen from the singular-value spectrum (values above 1e-3 of the largest). Capped to the pencil dimensions. pencil_factor : float Pencil parameter as a fraction of the sample count, L = round(factor·N); Hua & Sarkar recommend 0.33–0.5 for best noise rejection. Must lie in (0, 1) and yield 1 ≤ L ≤ N − 1. damping_threshold : float Damping ratio below which a mode is flagged poorly_damped. energy_floor : float Modes whose amplitude is below energy_floor times the largest modal amplitude are discarded as numerical noise.

Returns

tuple[OscillationMode, ...] The recovered modes, ordered by descending amplitude. Conjugate pairs are merged into one positive-frequency mode; near-zero-frequency real modes (pure decay) are reported with frequency_hz = 0.

Raises

ValueError If the signal, sampling rate, or parameters are invalid.

Source code in src/scpn_phase_orchestrator/monitor/oscillation_modes.py

def estimate_oscillation_modes(
    signal: FloatArray,
    fs: float,
    *,
    model_order: int | None = None,
    pencil_factor: float = 0.4,
    damping_threshold: float = DEFAULT_DAMPING_THRESHOLD,
    energy_floor: float = 1.0e-3,
) -> tuple[OscillationMode, ...]:
    """Estimate damped oscillation modes of a ringdown by the matrix-pencil method.

    Parameters
    ----------
    signal : FloatArray
        Real ringdown samples, uniformly sampled, length ``≥ 4``.
    fs : float
        Sampling frequency in hertz (``> 0``).
    model_order : int | None
        Number of modes (real exponentials plus conjugate pairs) to recover. When
        ``None`` the order is chosen from the singular-value spectrum (values
        above ``1e-3`` of the largest). Capped to the pencil dimensions.
    pencil_factor : float
        Pencil parameter as a fraction of the sample count, ``L = round(factor·N)``;
        Hua & Sarkar recommend ``0.33``–``0.5`` for best noise rejection. Must lie
        in ``(0, 1)`` and yield ``1 ≤ L ≤ N − 1``.
    damping_threshold : float
        Damping ratio below which a mode is flagged ``poorly_damped``.
    energy_floor : float
        Modes whose amplitude is below ``energy_floor`` times the largest modal
        amplitude are discarded as numerical noise.

    Returns
    -------
    tuple[OscillationMode, ...]
        The recovered modes, ordered by descending amplitude. Conjugate pairs are
        merged into one positive-frequency mode; near-zero-frequency real modes
        (pure decay) are reported with ``frequency_hz = 0``.

    Raises
    ------
    ValueError
        If the signal, sampling rate, or parameters are invalid.
    """
    samples = _validate_signal(signal)
    sample_rate = _positive_real(fs, "fs")
    threshold = _real_scalar(damping_threshold, "damping_threshold")
    floor = _non_negative_real(energy_floor, "energy_floor")
    n = samples.shape[0]
    pencil = _validated_pencil(pencil_factor, n)

    poles = _matrix_pencil_poles(samples, pencil, _validated_order(model_order))
    if poles.size == 0:
        return ()
    amplitudes, phases = _modal_residues(samples, poles)
    modes = _assemble_modes(poles, amplitudes, phases, sample_rate, threshold)
    return _prune_and_sort(modes, floor)