gds-domains/symbolic: Padé approximation for time delay modeling

[rororowyourboat]

Summary

Add delay.py to gds_domains/symbolic/ providing Padé approximation of pure time delays as rational transfer functions. Uses SymPy only.

Parent issue: #198 (classical control theory stack)
Depends on: #199 (transfer functions)

Motivation

Real control systems have transport delays, communication latency, and sampling delays. A pure delay e^{-sτ} cannot be represented as a finite-dimensional transfer function, but Padé approximation provides a rational polynomial approximation that preserves the phase characteristics. The ExecutionContract already has an observation_delay concept — this provides the frequency-domain model for analyzing its effects.

Proposed API

New file: gds_domains/symbolic/delay.py

def pade_approximation(delay: float, order: int = 3) -> "TransferFunction":
    """Padé approximation of e^{-sτ} as a rational transfer function.
    
    Order N produces an (N,N) approximation — N zeros, N poles, all-pass.
    Uses the closed-form Padé coefficients via SymPy.
    
    Args:
        delay: Time delay τ in seconds (must be > 0)
        order: Approximation order (1..8 typical; higher = more accurate
               but adds RHP zeros that complicate root locus)
    
    Returns:
        TransferFunction with num/den coefficients
    """

def delay_system(
    tf: "TransferFunction",
    delay: float,
    order: int = 3,
) -> "TransferFunction":
    """Cascade a transfer function with a Padé-approximated delay.
    
    Returns H(s) * Padé(s, τ) as a single TransferFunction by
    convolving the numerator/denominator polynomials.
    """

Implementation Notes

Padé formula

The (n, n) Padé approximant of e^{-sτ} is:

P_n(s) = Σ_{k=0}^{n} (-1)^k * C(n,k) * (sτ)^k / C(2n,k) * k!
Q_n(s) = Σ_{k=0}^{n} C(n,k) * (sτ)^k / C(2n,k) * k!

e^{-sτ} ≈ P_n(s) / Q_n(s)

SymPy has sympy.functions.combinatorial.numbers.nC and factorial for the coefficients. Alternatively, sympy.Rational for exact coefficient computation before converting to float.

Properties to validate

• All-pass: |H(jω)| = 1 for all ω (within numerical tolerance)
• Phase: matches e^{-jωτ} phase up to order 2N+1 in the Taylor expansion
• RHP zeros: order-N approximation has N zeros in the right half-plane (this is expected, not a bug)

Polynomial convolution

delay_system() multiplies two transfer functions by convolving coefficient lists:

result_num = numpy.convolve(tf.num, pade.num)  # or pure-Python convolution
result_den = numpy.convolve(tf.den, pade.den)

Since this is in the symbolic layer (no numpy), use a pure-Python polynomial multiplication or SymPy's Poly.mul().

Key Files

• packages/gds-domains/gds_domains/symbolic/transfer.py — TransferFunction (from gds-domains/symbolic: transfer functions, poles/zeros, and sensitivity functions #199)
• packages/gds-domains/gds_domains/symbolic/linearize.py — reference for data conventions
• packages/gds-framework/gds/execution.py — ExecutionContract.observation_delay (conceptual link)

Testing

• test_delay.py:
  • Order-1 Padé of 0.1s delay: verify |H(jω)| ≈ 1 at 10 frequencies
  • Order-3 Padé: verify phase matches -ωτ to within 1% up to ω = 5/τ
  • delay_system() with 1/(s+1) and 0.5s delay: verify pole count = original + order
  • Edge case: delay=0 returns identity TF 1/1

Concepts Addressed (MATLAB Tech Talks)

• Video 11: Time delay — transport delay, latency, stability margin effects
• Video 12: Padé approximations — rational polynomial delay models, all-pass filters
• Video 13: Non-minimum phase — Padé introduces RHP zeros (expected)