Modal Participation — Damping Controllability and Mode Shapes

monitor.modal_participation is the model-based companion to the data-driven oscillation-mode estimator. The matrix-pencil estimator reads a measured ringdown and returns each mode's damping; this module reads the network model and returns the two further things a damping-recommendation engine needs: where an inter-area mode swings (mode shape and participation factors) and which actuator damps it best (modal controllability). A single ringdown signal cannot give those — they come from the eigenstructure of the small-signal state matrix.

The small-signal state matrix

phase_network_jacobian linearises the engine's Sakaguchi–Kuramoto coupling dynamics about an operating point θ* (typically a synchronised fixed point reached by running the engine):

• off-diagonal: J_ik = K_ik cos(θ_k* − θ_i* − α_ik)
• diagonal: J_ii = −Σ_{k≠i} K_ik cos(θ_k* − θ_i* − α_ik) − ζ cos(Ψ − θ_i*)

This matches the integrator's derivative exactly. For a symmetric, lag-free network the Jacobian is symmetric negative-semidefinite — overdamped, with real eigenvalues and no oscillation. The oscillatory inter-area modes that reliability standards screen for appear precisely when the Sakaguchi phase lag α or a directed (asymmetric) coupling break that symmetry, which is the regime the SCPN engine runs in. A second-order swing/companion state matrix can be analysed directly too — analyse_network_modes accepts any continuous-time matrix.

Modal analysis

analyse_network_modes eigen-decomposes the state matrix A (LAPACK geev via NumPy) and, for each eigenvalue λ = σ + jω, reports:

1. Frequency f = |ω| / 2π and damping ratio ζ = −σ / |λ| (Kundur 1994, §12). A growing (unstable) mode has a negative damping ratio; the marginal global-phase mode of a Kuramoto Jacobian appears at f = 0, ζ = 0.
2. Mode shape — the right eigenvector φ_i, unit-norm and phase-anchored, so |φ_ki| is oscillator k's relative swing amplitude and ∠φ_ki its relative phase. Anti-phase entries identify the machines that swing against each other.
3. Participation factors p_ki = φ_ki · ψ_ik (Pérez-Arriaga, Verghese & Schweppe 1982), with ψ_i the matching left eigenvector — dimensionless, real, non-negative, summing to one. participation[k] is how much oscillator k shapes the mode; dominant_state is the largest.
4. Modal controllability |ψ_i · b_j| for each column b_j of an optional input matrix B — how strongly input j can drive the mode. dominant_input is the most effective actuator: where to add damping.

Complex eigenvalues of a real matrix occur in conjugate pairs; only the non-negative-frequency member of each pair is reported. Modes are ordered by ascending damping ratio — the least-damped, most critical mode first. A defective (non-diagonalisable) state matrix is rejected: its left eigenvectors, and the participation factors built from them, do not exist.

All public state, input, coupling, phase, and phase-lag arrays must contain finite real numeric evidence before conversion. Boolean, complex, numeric-string, overflow, and broken array-protocol payloads fail closed; real numeric object arrays remain valid. A directly constructed NetworkMode independently replays eigenvalue-derived frequency and damping, canonical unit/phase anchoring of its mode shape, normalised participation and dominant state, controllability/dominant-input pairing, and primitive field types. It owns read-only copies of all array evidence, so later caller mutation cannot alter a published modal record.

Relationship to oscillation_modes

The two pieces describe the same linear modes from opposite ends. For a linear ringdown x(t) = e^{At} x₀, estimate_oscillation_modes recovers the eigenvalues of A from the measured signal, and analyse_network_modes(A) computes them from the model — they agree (this cross-check is part of the test suite). Use the estimator when you only have measurements; use this module when you have the model and need mode shape, participation, or controllability.

NumPy floor

The analysis is one offline eigen-decomposition of a modest state matrix (LAPACK via NumPy), not a per-step hot path, so it stays on the NumPy floor — the same judgement as oscillation_modes and autotune.freq_id; it carries no multi-language acceleration chain.

Review-only

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

modal_participation

Model-based modal participation and damping controllability of the phase network.

Where :mod:~scpn_phase_orchestrator.monitor.oscillation_modes recovers modal damping from a measured ringdown (data-driven), this module answers the two questions a damping-recommendation engine needs from the network model: for a poorly-damped inter-area mode, which oscillators swing in it (mode shape and participation factors → where the mode is observable) and which actuators can damp it (modal controllability → where to act). A single ringdown signal cannot give those — they come from the eigenstructure of the small-signal state matrix.

The small-signal state matrix of the Sakaguchi–Kuramoto network is the Jacobian of the engine's coupling dynamics θ̇_i = ω_i + Σ_j K_ij sin(θ_j − θ_i − α_ij) + ζ sin(Ψ − θ_i) linearised about an operating point θ* (:func:phase_network_jacobian): J_ik = K_ik cos(θ_k* − θ_i* − α_ik) off the diagonal and J_ii = −Σ_{k≠i} K_ik cos(θ_k* − θ_i* − α_ik) − ζ cos(Ψ − θ_i*) on it. For a symmetric, lag-free network the Jacobian is symmetric negative-semidefinite (overdamped, real eigenvalues); the oscillatory inter-area modes that reliability standards screen for appear precisely when the phase lag α or a directed (asymmetric) coupling break that symmetry — the regime the SCPN engine runs in.

:func:analyse_network_modes performs the small-signal modal analysis of any continuous-time state matrix A (so a second-order swing/companion form can be fed in directly too). It eigen-decomposes A (LAPACK geev via NumPy), reads each eigenvalue λ = σ + jω as a mode of frequency f = |ω| / 2π and damping ratio ζ = −σ / |λ| (Kundur 1994, §12), and from the right eigenvector φ_i and the matching left eigenvector ψ_i (the rows of A's eigenvector inverse) forms the dimensionless participation factor p_ki = φ_ki · ψ_ik (Pérez-Arriaga, Verghese & Schweppe 1982) and the modal controllability |ψ_i · b_j| of mode i from input j. Complex eigenvalues of a real matrix occur in conjugate pairs; only the non-negative-frequency member of each pair is reported. The analysis is diagnostic only — it reads a model and reports modes; it never changes bindings, layers, or coupling.

The whole analysis is one offline eigen-decomposition of a modest state matrix (LAPACK via NumPy), not a per-step hot path, so it stays on the NumPy floor — the same judgement as oscillation_modes and autotune.freq_id.

References

• Kundur, P. 1994, Power System Stability and Control (McGraw-Hill), §12 — small-signal stability, eigenvalues, mode shapes, participation factors.
• Pérez-Arriaga, I. J., Verghese, G. C. & Schweppe, F. C. 1982, IEEE Trans. Power App. Syst. PAS-101(9):3117–3125 — selective modal analysis; participation factors.
• Dörfler, F. & Bullo, F. 2014, Automatica 50(6):1539–1564 — synchronisation in networks of phase oscillators; the Kuramoto stability Jacobian.
• NERC PRC-028 (oscillation monitoring) — damping-ratio screening of inter-area modes.

Classes

NetworkMode dataclass

NetworkMode(
    eigenvalue: complex,
    frequency_hz: float,
    damping_ratio: float,
    mode_shape: ComplexArray,
    participation: FloatArray,
    dominant_state: int,
    controllability: FloatArray | None,
    dominant_input: int | None,
    poorly_damped: bool,
)

One small-signal mode of a phase-oscillator network.

Equality is identity-based (eq=False): the array fields make value equality ambiguous, and modes are compared field by field, never as wholes.

Attributes

eigenvalue : complex The continuous-time eigenvalue σ + jω (rad/s) of the state matrix, taken from the non-negative-frequency member of its conjugate pair. frequency_hz : float Modal oscillation frequency |ω| / 2π in hertz (≥ 0); 0 for a non-oscillatory (real-eigenvalue) mode. damping_ratio : float Dimensionless damping ratio ζ = −σ / |λ|; > 0 is stable, < 0 is growing (unstable), 0 is the marginal global-phase mode. mode_shape : ComplexArray Right eigenvector, unit Euclidean norm and phase-anchored so the largest-magnitude entry is real-positive: per-oscillator relative amplitude (|·|) and relative phase (∠) of the swing. participation : FloatArray Real participation factors over oscillators, ≥ 0 and summing to 1; participation[k] measures how much oscillator k shapes the mode. dominant_state : int Index of the oscillator with the largest participation factor. controllability : FloatArray | None Per-input modal controllability |ψ_i^{H} b_j| when an input matrix is supplied, otherwise None; larger means input j damps the mode more. dominant_input : int | None Index of the most effective input, or None without an input matrix. poorly_damped : bool Whether damping_ratio is below the screening threshold.

Methods:

__post_init__

__post_init__() -> None

Validate, normalise, and own the published modal evidence.

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

def __post_init__(self) -> None:
    """Validate, normalise, and own the published modal evidence."""
    eigenvalue = _complex_scalar(self.eigenvalue, "eigenvalue")
    frequency = _nonnegative_scalar(self.frequency_hz, "frequency_hz")
    damping = _real_scalar(self.damping_ratio, "damping_ratio")
    expected_frequency = abs(eigenvalue.imag) / (2.0 * np.pi)
    magnitude = abs(eigenvalue)
    expected_damping = (
        -eigenvalue.real / magnitude if magnitude > _ZERO_EIGENVALUE else 0.0
    )
    if not np.isclose(
        frequency,
        expected_frequency,
        rtol=8.0 * np.finfo(np.float64).eps,
        atol=8.0 * np.finfo(np.float64).eps,
    ):
        raise ValueError("frequency_hz must match the eigenvalue")
    if not np.isclose(
        damping,
        expected_damping,
        rtol=8.0 * np.finfo(np.float64).eps,
        atol=8.0 * np.finfo(np.float64).eps,
    ):
        raise ValueError("damping_ratio must match the eigenvalue")

    mode_shape = _validate_complex_vector(self.mode_shape, "mode_shape")
    norm = float(np.linalg.norm(mode_shape))
    if not np.isclose(norm, 1.0, rtol=0.0, atol=1.0e-12):
        raise ValueError("mode_shape must have unit Euclidean norm")
    anchor = mode_shape[int(np.argmax(np.abs(mode_shape)))]
    if abs(float(anchor.imag)) > 1.0e-12 or float(anchor.real) <= 0.0:
        raise ValueError("mode_shape must be anchored real-positive")

    participation = _validate_real_array(self.participation, "participation")
    if participation.ndim != 1 or participation.shape != mode_shape.shape:
        raise ValueError("participation must match the one-dimensional mode_shape")
    if np.any(participation < 0.0) or not np.isclose(
        float(participation.sum()),
        1.0,
        rtol=0.0,
        atol=1.0e-12,
    ):
        raise ValueError("participation must be non-negative and sum to one")
    dominant_state = _validate_index(
        self.dominant_state,
        "dominant_state",
        participation.size,
    )
    if dominant_state != int(np.argmax(participation)):
        raise ValueError("dominant_state must identify maximum participation")

    controllability: FloatArray | None
    dominant_input: int | None
    if self.controllability is None:
        controllability = None
        if self.dominant_input is not None:
            raise ValueError("dominant_input must be None without controllability")
        dominant_input = None
    else:
        controllability = _validate_real_array(
            self.controllability,
            "controllability",
        )
        if controllability.ndim != 1 or controllability.size == 0:
            raise ValueError("controllability must be a non-empty vector")
        if np.any(controllability < 0.0):
            raise ValueError("controllability must be non-negative")
        if self.dominant_input is None:
            raise ValueError("dominant_input is required with controllability")
        dominant_input = _validate_index(
            self.dominant_input,
            "dominant_input",
            controllability.size,
        )
        if dominant_input != int(np.argmax(controllability)):
            raise ValueError("dominant_input must identify maximum controllability")
    poorly_damped = _plain_bool(self.poorly_damped, "poorly_damped")
    mode_shape.setflags(write=False)
    participation.setflags(write=False)
    if controllability is not None:
        controllability.setflags(write=False)
    object.__setattr__(self, "eigenvalue", eigenvalue)
    object.__setattr__(self, "frequency_hz", frequency)
    object.__setattr__(self, "damping_ratio", damping)
    object.__setattr__(self, "mode_shape", mode_shape)
    object.__setattr__(self, "participation", participation)
    object.__setattr__(self, "dominant_state", dominant_state)
    object.__setattr__(self, "controllability", controllability)
    object.__setattr__(self, "dominant_input", dominant_input)
    object.__setattr__(self, "poorly_damped", poorly_damped)

to_dict

to_dict() -> dict[str, object]

Return a JSON-serialisable mapping of the mode.

Returns

dict[str, object] The eigenvalue (as [real, imag]), frequency, damping ratio, mode shape (as [real, imag] pairs), participation factors, dominant oscillator, per-input controllability (or None), dominant input, and the poorly-damped flag.

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

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

    Returns
    -------
    dict[str, object]
        The eigenvalue (as ``[real, imag]``), frequency, damping ratio, mode
        shape (as ``[real, imag]`` pairs), participation factors, dominant
        oscillator, per-input controllability (or ``None``), dominant input,
        and the poorly-damped flag.
    """
    controllability = (
        None
        if self.controllability is None
        else [float(value) for value in self.controllability]
    )
    return {
        "eigenvalue": [float(self.eigenvalue.real), float(self.eigenvalue.imag)],
        "frequency_hz": self.frequency_hz,
        "damping_ratio": self.damping_ratio,
        "mode_shape": [
            [float(value.real), float(value.imag)] for value in self.mode_shape
        ],
        "participation": [float(value) for value in self.participation],
        "dominant_state": self.dominant_state,
        "controllability": controllability,
        "dominant_input": self.dominant_input,
        "poorly_damped": self.poorly_damped,
    }

Functions:

analyse_network_modes

analyse_network_modes(
    state_matrix: FloatArray,
    *,
    input_matrix: FloatArray | None = None,
    damping_threshold: float = DEFAULT_DAMPING_THRESHOLD,
) -> tuple[NetworkMode, ...]

Decompose a continuous-time state matrix into damped modes with participation.

Parameters

state_matrix : FloatArray Real square Jacobian A = ∂ẋ/∂x of the linearised dynamics, shape (N, N); build it from a phase network with :func:phase_network_jacobian. input_matrix : FloatArray | None Real input matrix B, shape (N, M), mapping M actuator inputs into the state derivative; when given, each mode reports per-input modal controllability. None skips the controllability analysis. damping_threshold : float Damping ratio below which a mode is flagged poorly_damped.

Returns

tuple[NetworkMode, ...] One mode per non-negative-frequency eigenvalue, ordered by ascending damping ratio (least-damped, most critical first). The marginal global-phase mode of a Kuramoto Jacobian appears with frequency_hz = 0 and damping_ratio = 0.

Raises

ValueError If the state matrix or input matrix is invalid, or the state matrix is defective (not diagonalisable, so left eigenvectors do not exist).

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

def analyse_network_modes(
    state_matrix: FloatArray,
    *,
    input_matrix: FloatArray | None = None,
    damping_threshold: float = DEFAULT_DAMPING_THRESHOLD,
) -> tuple[NetworkMode, ...]:
    """Decompose a continuous-time state matrix into damped modes with participation.

    Parameters
    ----------
    state_matrix : FloatArray
        Real square Jacobian ``A = ∂ẋ/∂x`` of the linearised dynamics, shape
        ``(N, N)``; build it from a phase network with
        :func:`phase_network_jacobian`.
    input_matrix : FloatArray | None
        Real input matrix ``B``, shape ``(N, M)``, mapping ``M`` actuator inputs
        into the state derivative; when given, each mode reports per-input modal
        controllability. ``None`` skips the controllability analysis.
    damping_threshold : float
        Damping ratio below which a mode is flagged ``poorly_damped``.

    Returns
    -------
    tuple[NetworkMode, ...]
        One mode per non-negative-frequency eigenvalue, ordered by ascending
        damping ratio (least-damped, most critical first). The marginal
        global-phase mode of a Kuramoto Jacobian appears with
        ``frequency_hz = 0`` and ``damping_ratio = 0``.

    Raises
    ------
    ValueError
        If the state matrix or input matrix is invalid, or the state matrix is
        defective (not diagonalisable, so left eigenvectors do not exist).
    """
    matrix = _validate_square_matrix(state_matrix, "state_matrix")
    threshold = _real_scalar(damping_threshold, "damping_threshold")
    inputs = _validate_input_matrix(input_matrix, matrix.shape[0])

    raw_eigenvalues, raw_right = np.linalg.eig(matrix)
    eigenvalues = raw_eigenvalues.astype(np.complex128)
    right = raw_right.astype(np.complex128)
    # A defective (non-diagonalisable) matrix has a rank-deficient eigenvector
    # matrix, so its left eigenvectors — and the participation factors built from
    # them — do not exist.
    if np.linalg.matrix_rank(right) < right.shape[0]:
        raise ValueError("state_matrix is defective (not diagonalisable)")
    left = np.linalg.inv(right).astype(np.complex128)

    modes = [
        _build_mode(index, eigenvalues, right, left, inputs, threshold)
        for index in range(eigenvalues.shape[0])
        if eigenvalues[index].imag >= 0.0
    ]
    modes.sort(key=lambda mode: mode.damping_ratio)
    return tuple(modes)

phase_network_jacobian

phase_network_jacobian(
    coupling: FloatArray,
    phases: FloatArray,
    *,
    phase_lag: FloatArray | None = None,
    drive_strength: float = 0.0,
    drive_phase: float = 0.0,
) -> FloatArray

Build the Sakaguchi–Kuramoto small-signal Jacobian at an operating point.

The Jacobian is the state matrix of the engine's coupling dynamics linearised about phases; feed it to :func:analyse_network_modes. It matches the integrator's derivative exactly: J_ik = K_ik cos(θ_k − θ_i − α_ik) off the diagonal and J_ii = −Σ_{k≠i} K_ik cos(θ_k − θ_i − α_ik) − ζ cos(Ψ − θ_i).

Parameters

coupling : FloatArray Coupling matrix K_nm, shape (N, N), with a zero diagonal (no self-coupling), as the integrator requires. phases : FloatArray Operating-point phases θ* in radians, shape (N,) — typically a synchronised fixed point reached by running the engine. phase_lag : FloatArray | None Sakaguchi phase-lag matrix α in radians, shape (N, N); None means no lag (zeros). drive_strength : float External-drive strength ζ; the default 0 gives the free-network Jacobian whose modes are the inter-area oscillations. drive_phase : float External-drive reference phase Ψ in radians (used only when drive_strength is non-zero).

Returns

FloatArray The (N, N) small-signal Jacobian.

Raises

ValueError If the coupling matrix, phases, or phase-lag matrix are invalid, or the coupling diagonal is non-zero.

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

def phase_network_jacobian(
    coupling: FloatArray,
    phases: FloatArray,
    *,
    phase_lag: FloatArray | None = None,
    drive_strength: float = 0.0,
    drive_phase: float = 0.0,
) -> FloatArray:
    """Build the Sakaguchi–Kuramoto small-signal Jacobian at an operating point.

    The Jacobian is the state matrix of the engine's coupling dynamics linearised
    about ``phases``; feed it to :func:`analyse_network_modes`. It matches the
    integrator's derivative exactly: ``J_ik = K_ik cos(θ_k − θ_i − α_ik)`` off the
    diagonal and ``J_ii = −Σ_{k≠i} K_ik cos(θ_k − θ_i − α_ik) − ζ cos(Ψ − θ_i)``.

    Parameters
    ----------
    coupling : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``, with a zero diagonal (no
        self-coupling), as the integrator requires.
    phases : FloatArray
        Operating-point phases ``θ*`` in radians, shape ``(N,)`` — typically a
        synchronised fixed point reached by running the engine.
    phase_lag : FloatArray | None
        Sakaguchi phase-lag matrix ``α`` in radians, shape ``(N, N)``; ``None``
        means no lag (zeros).
    drive_strength : float
        External-drive strength ``ζ``; the default ``0`` gives the free-network
        Jacobian whose modes are the inter-area oscillations.
    drive_phase : float
        External-drive reference phase ``Ψ`` in radians (used only when
        ``drive_strength`` is non-zero).

    Returns
    -------
    FloatArray
        The ``(N, N)`` small-signal Jacobian.

    Raises
    ------
    ValueError
        If the coupling matrix, phases, or phase-lag matrix are invalid, or the
        coupling diagonal is non-zero.
    """
    matrix = _validate_square_matrix(coupling, "coupling")
    n = matrix.shape[0]
    angles = _validate_vector(phases, "phases", n)
    lag = _validate_phase_lag(phase_lag, n)
    strength = _real_scalar(drive_strength, "drive_strength")
    reference = _real_scalar(drive_phase, "drive_phase")
    if not np.allclose(np.diag(matrix), 0.0, rtol=0.0, atol=1.0e-15):
        raise ValueError("coupling self-coupling diagonal must be zero")

    differences = angles[np.newaxis, :] - angles[:, np.newaxis] - lag
    jacobian = matrix * np.cos(differences)
    row_sums = jacobian.sum(axis=1)
    drive = strength * np.cos(reference - angles)
    np.fill_diagonal(jacobian, -(row_sums + drive))
    return np.ascontiguousarray(jacobian, dtype=np.float64)