gds-domains/symbolic: transfer functions, poles/zeros, and sensitivity functions

[rororowyourboat]

Summary

Add transfer.py to gds_domains/symbolic/ providing symbolic transfer function representation, pole/zero computation, controllability/observability matrices, non-minimum phase detection, and sensitivity functions (Gang of Six). Uses SymPy only — no scipy or numpy.

Parent issue: #198 (classical control theory stack)

Motivation

LinearizedSystem already provides (A, B, C, D) as list[list[float]]. The transfer function H(s) = C(sI - A)^{-1}B + D is a one-hop derivation that unlocks the entire frequency-domain analysis chain. Without it, users have to leave the GDS ecosystem for basic control analysis.

Proposed API

New file: gds_domains/symbolic/transfer.py

from dataclasses import dataclass

@dataclass(frozen=True)
class TransferFunction:
    """SISO transfer function as coefficient lists (descending powers of s)."""
    num: list[float]       # numerator coefficients [b_n, b_{n-1}, ..., b_0]
    den: list[float]       # denominator coefficients [a_m, a_{m-1}, ..., a_0]
    input_name: str = ""
    output_name: str = ""

@dataclass(frozen=True)
class TransferFunctionMatrix:
    """MIMO transfer function matrix (p outputs × m inputs)."""
    elements: list[list[TransferFunction]]
    input_names: list[str]
    output_names: list[str]

def ss_to_tf(ls: "LinearizedSystem") -> TransferFunctionMatrix:
    """Convert state-space (A, B, C, D) to transfer function matrix.
    
    H(s) = C(sI - A)^{-1}B + D
    Uses SymPy's Matrix.adjugate() and det() for symbolic computation,
    then extracts polynomial coefficients.
    """

def characteristic_polynomial(ls: "LinearizedSystem") -> list[float]:
    """Compute det(sI - A) as coefficient list."""

def poles(tf: TransferFunction) -> list[complex]:
    """Roots of denominator polynomial."""

def zeros(tf: TransferFunction) -> list[complex]:
    """Roots of numerator polynomial."""

def is_minimum_phase(tf: TransferFunction) -> bool:
    """True if all zeros have Re(z) < 0."""

def controllability_matrix(ls: "LinearizedSystem") -> list[list[float]]:
    """C_c = [B, AB, A^2B, ..., A^{n-1}B]"""

def observability_matrix(ls: "LinearizedSystem") -> list[list[float]]:
    """C_o = [C; CA; CA^2; ...; CA^{n-1}]"""

def is_controllable(ls: "LinearizedSystem") -> bool:
    """rank(C_c) == n"""

def is_observable(ls: "LinearizedSystem") -> bool:
    """rank(C_o) == n"""

def sensitivity(plant: TransferFunction, controller: TransferFunction) -> dict[str, TransferFunction]:
    """Gang of Six sensitivity functions.
    
    Returns dict with keys: S, T, CS, PS, KS, KPS
    - S  = 1 / (1 + L)           — sensitivity
    - T  = L / (1 + L)           — complementary sensitivity
    - CS = C / (1 + L)           — control sensitivity
    - PS = P / (1 + L)           — load sensitivity
    - KS = K / (1 + L)           — noise sensitivity (K=controller)
    - KPS = KP / (1 + L)         — input disturbance sensitivity
    where L = P * K (loop transfer function)
    """

Export from __init__.py

Add TransferFunction, TransferFunctionMatrix, ss_to_tf, poles, zeros to gds_domains.symbolic public API.

Implementation Notes

SymPy approach

import sympy
s = sympy.Symbol("s")
A_mat = sympy.Matrix(ls.A)
I_n = sympy.eye(len(ls.A))
# H(s) = C @ (s*I - A).inv() @ B + D
# Use adjugate/det to avoid symbolic inversion
det_sIA = (s * I_n - A_mat).det()
adj_sIA = (s * I_n - A_mat).adjugate()
# Each element H_ij(s) = (C[i,:] @ adj @ B[:,j]) / det + D[i,j]

This avoids numerical issues with symbolic matrix inversion for larger systems.

Interface contract

TransferFunction uses list[float] for coefficients, matching the list[list[float]] convention used by LinearizedSystem. No numpy arrays in the public API.

Functions that accept LinearizedSystem do so by type annotation — the actual access is to .A, .B, .C, .D fields, so any object with those list[list[float]] attributes would work.

Dependency

SymPy only — already the [symbolic] extra. No new dependencies.

Key Files

• packages/gds-domains/gds_domains/symbolic/linearize.py — LinearizedSystem definition (input to this work)
• packages/gds-domains/gds_domains/symbolic/model.py — SymbolicControlModel.linearize() dispatch
• packages/gds-domains/gds_domains/symbolic/__init__.py — public API exports
• packages/gds-domains/gds_domains/symbolic/hamiltonian.py — reference for SymPy usage patterns in this package

Testing

• test_transfer.py with known analytic cases:
  • 1st-order system: H(s) = 1/(s+1) from A=[-1], B=[1], C=[1], D=[0]
  • 2nd-order underdamped: known poles at −ζω_n ± jω_n√(1−ζ²)
  • Double integrator: H(s) = 1/s² with two poles at origin
  • Controllability/observability for observable canonical form (full rank) vs. uncontrollable system (rank-deficient)
  • Non-minimum phase system with known RHP zero
  • Gang of Six: verify S + T = 1 identity

Concepts Addressed (MATLAB Tech Talks)

• Video 2: Transfer functions, poles/zeros, S-domain
• Video 4: System stability via poles (precursor to Bode/Nyquist)
• Video 7: Gang of Six sensitivity functions
• Video 13: Non-minimum phase detection (RHP zeros)