gds-analysis: numerical linear systems analysis, discretization, and controller synthesis

[rororowyourboat]

Summary

Add linear.py to gds_analysis/ providing numerical linear systems analysis: eigenvalue stability, frequency response computation, stability margins, discretization (c2d), LQR synthesis, Kalman filter design, and gain scheduling. All behind the existing [continuous] extra (scipy + numpy already available).

Parent issue: #198 (classical control theory stack)
Depends on: #199 (transfer functions — for margin computation from TF coefficients)

Motivation

gds-domains/symbolic computes LinearizedSystem (A, B, C, D) and transfer functions symbolically. But users need fast numerical answers: "is this stable?" (eigenvalues), "what are the margins?" (frequency response), "what gains should I use?" (LQR), "how do I implement this digitally?" (discretization). These are scipy one-liners on the matrices we already have — the analysis package is the natural home.

Proposed API

New file: gds_analysis/linear.py (behind [continuous] extra)

All functions accept list[list[float]] matrices — no LinearizedSystem import, no coupling to gds-domains.

Stability Analysis

def eigenvalues(A: list[list[float]]) -> list[complex]:
    """Eigenvalues of state matrix A. Wraps numpy.linalg.eigvals."""

def is_stable(A: list[list[float]], continuous: bool = True) -> bool:
    """Check asymptotic stability.
    
    Continuous: all Re(λ) < 0
    Discrete:   all |λ| < 1
    """

def is_marginally_stable(A: list[list[float]], continuous: bool = True) -> bool:
    """Check marginal stability (eigenvalues on boundary, none in RHP/outside unit circle)."""

Frequency Response

def frequency_response(
    A: list[list[float]],
    B: list[list[float]],
    C: list[list[float]],
    D: list[list[float]],
    omega: list[float] | None = None,
) -> tuple[list[float], list[float], list[float]]:
    """Compute frequency response H(jω) numerically.
    
    Returns (omega, magnitude_db, phase_deg) as plain lists.
    If omega is None, auto-selects frequency range from eigenvalues.
    Wraps scipy.signal.bode or scipy.signal.freqresp.
    """

def gain_margin(num: list[float], den: list[float]) -> tuple[float, float]:
    """Gain margin (dB) and frequency (rad/s) where phase = -180°.
    Wraps scipy.signal.margin or custom phase-crossing detection.
    """

def phase_margin(num: list[float], den: list[float]) -> tuple[float, float]:
    """Phase margin (deg) and gain crossover frequency (rad/s).
    Wraps scipy.signal.margin.
    """

Discretization

def discretize(
    A: list[list[float]],
    B: list[list[float]],
    C: list[list[float]],
    D: list[list[float]],
    dt: float,
    method: str = "tustin",
) -> tuple[list[list[float]], list[list[float]], list[list[float]], list[list[float]]]:
    """Convert continuous state-space to discrete: (A, B, C, D) → (Ad, Bd, Cd, Dd).
    
    Methods: "tustin" (bilinear), "zoh" (zero-order hold), 
             "euler" (forward Euler), "backward_euler"
    
    Wraps scipy.signal.cont2discrete.
    Returns (Ad, Bd, Cd, Dd) as list[list[float]].
    """

Controller Synthesis

def lqr(
    A: list[list[float]],
    B: list[list[float]],
    Q: list[list[float]],
    R: list[list[float]],
    N: list[list[float]] | None = None,
) -> tuple[list[list[float]], list[list[float]], list[complex]]:
    """Linear Quadratic Regulator.
    
    Solves the continuous algebraic Riccati equation:
        A'P + PA - PBR^{-1}B'P + Q = 0
    
    Returns (K, P, E):
        K — optimal gain matrix (u = -Kx)
        P — Riccati solution
        E — closed-loop eigenvalues
    
    Wraps scipy.linalg.solve_continuous_are.
    """

def dlqr(
    Ad: list[list[float]],
    Bd: list[list[float]],
    Q: list[list[float]],
    R: list[list[float]],
) -> tuple[list[list[float]], list[list[float]], list[complex]]:
    """Discrete-time LQR. Wraps scipy.linalg.solve_discrete_are."""

def kalman(
    A: list[list[float]],
    C: list[list[float]],
    Q_process: list[list[float]],
    R_measurement: list[list[float]],
) -> tuple[list[list[float]], list[list[float]]]:
    """Steady-state Kalman filter (observer) gain.
    
    Dual of LQR: solves ARE for the observer.
    Returns (L, P): observer gain and error covariance.
    """

def gain_schedule(
    linearize_fn: "Callable[[dict], tuple[list[list[float]], ...]]",
    operating_points: list[dict],
    Q: list[list[float]],
    R: list[list[float]],
) -> list[tuple[dict, list[list[float]], list[complex]]]:
    """Compute LQR gains at multiple operating points.
    
    linearize_fn: takes operating point dict, returns (A, B, C, D)
    Returns list of (point, K, closed_loop_eigenvalues) tuples.
    """

Implementation Notes

scipy functions used

GDS function	scipy wrapper	Module
eigenvalues	numpy.linalg.eigvals	numpy
frequency_response	scipy.signal.bode / scipy.signal.freqresp	scipy.signal
gain_margin / phase_margin	scipy.signal.margin or custom	scipy.signal
discretize	scipy.signal.cont2discrete	scipy.signal
lqr	scipy.linalg.solve_continuous_are	scipy.linalg
dlqr	scipy.linalg.solve_discrete_are	scipy.linalg
kalman	dual of solve_continuous_are	scipy.linalg

Interface contract

• Input: list[list[float]] matrices, matching LinearizedSystem field types
• Output: list[list[float]], list[float], list[complex] — never numpy arrays
• Internal: convert to numpy at entry, compute, convert back at exit
• No gds-domains import: functions accept raw matrices, not typed objects

Dependency

All behind the existing [continuous] extra which already pulls gds-continuous[scipy]>=0.1.0 + numpy>=1.26. Consider adding a [control] extra alias for discoverability.

Error handling

• is_stable() on an empty matrix returns True (vacuously stable)
• lqr() raises ValueError if (A, B) is not stabilizable
• kalman() raises ValueError if (A, C) is not detectable
• discretize() raises ValueError for unknown method strings

Key Files

• packages/gds-analysis/gds_analysis/backward_reachability.py — reference for scipy usage patterns in this package
• packages/gds-analysis/gds_analysis/adapter.py — the discrete-time bridge (this is the continuous-time analog)
• packages/gds-domains/gds_domains/symbolic/linearize.py — LinearizedSystem definition (data source)
• packages/gds-continuous/gds_continuous/engine.py — scipy usage patterns

Testing

• test_linear.py:
  • Eigenvalue: 2×2 stable system A = [[-1, 0], [0, -2]] → eigenvalues {-1, -2}, is_stable=True
  • Eigenvalue: unstable system with positive eigenvalue → is_stable=False
  • Frequency response: 1/(s+1) → verify -3dB at ω=1, -45° at ω=1
  • Margins: known system with 6dB gain margin, 45° phase margin
  • Discretize: A=[-1], B=[1], C=[1], D=[0] with Tustin at dt=0.1 → verify against scipy.signal.cont2discrete directly
  • LQR: double integrator with identity Q, R → verify closed-loop eigenvalues in LHP
  • Kalman: dual of LQR test case
  • Gain scheduling: linearize spring-mass at 3 operating points → verify K varies

Concepts Addressed (MATLAB Tech Talks)

• Video 1: Stability analysis, LQR, Kalman filter, state estimation
• Video 2: Poles (eigenvalues), system stability
• Video 4: Frequency response, gain/phase margin
• Video 5: Discretization (Tustin, ZOH, Euler), state-space conversion
• Video 7: Stability margins for loop shaping
• Video 9: Gain scheduling via multi-point linearization
• Video 11: Stability margin effects of time delay
• Video 14: Z-transform (discretization bridge to discrete-time)

[Whatsonyourmind]

The proposed linear.py API is a clean fit for the symbolic LinearizedSystem output, and the "scipy one-liners on matrices we already have" framing is exactly right — most of the technical work is already done upstream. A few things worth thinking about before the first implementation, because the numerical quality of these operations depends more on preprocessing than on the scipy call itself.

1. LQR and Kalman filter design are mathematical duals — expose that structure in the API.

The continuous-time LQR problem solves the algebraic Riccati equation (CARE):

$$ A^T P + P A - P B R^{-1} B^T P + Q = 0 $$

for $P$ positive semi-definite, then constructs $K = R^{-1} B^T P$. The Kalman filter steady-state gain problem solves a mathematically identical CARE with substitutions:

$$ A^T \to A, \quad B \to C^T, \quad Q \to Q_w \text{ (process noise)}, \quad R \to R_v \text{ (measurement noise)} $$

and then constructs $L = P C^T R_v^{-1}$. scipy's solve_continuous_are(A, B, Q, R) handles both cases — for LQR you pass the arguments as-is, for Kalman design you pass (A.T, C.T, Q_w, R_v) and transpose the result. This means you can have a single internal _care_solve(...) helper and two public API functions that differ only in the argument shuffle.

A cleaner API would expose this explicitly:

def lqr(A, B, Q, R, N=None):
    """Continuous LQR. Returns (K, P, E).
    K = R^{-1} B^T P, u = -K x
    """
    P = _care_solve(A, B, Q, R, N)
    K = _solve_posdef(R, B.T @ P)  # more stable than R^{-1} @ B^T @ P
    E = eigenvalues(A - B @ K)
    return _as_list(K), _as_list(P), E

def kalman_design(A, C, Q_w, R_v):
    """Continuous steady-state Kalman filter gain. Returns (L, P, E).
    L = P C^T R_v^{-1}, x_hat' = A x_hat + L (y - C x_hat)
    """
    P = _care_solve(A.T, C.T, Q_w, R_v)
    L = _solve_posdef(R_v, C @ P)  # returns L^T; transpose
    L = L.T
    E = eigenvalues(A - L @ C)
    return _as_list(L), _as_list(P), E

Sharing the _care_solve helper means bugs in the underlying solve get fixed once, not twice, and users who understand the duality can pass control-synthesis-shaped inputs to kalman_design or vice versa without being surprised.

2. The LQR Riccati solve is ill-conditioned for slow-mode systems — use the Schur method explicitly.

scipy.linalg.solve_continuous_are has two solver paths: the default Schur decomposition method, and the older matrix-sign function method. For well-conditioned problems both work fine, but for systems with eigenvalues close to the imaginary axis (slow modes, common in cadCAD-style dynamical system specifications where mass-conservation implies a near-zero eigenvalue), the matrix-sign method can fail to converge or return a non-symmetric $P$.

The fix is to explicitly request the Schur method and verify the output is symmetric positive-definite before returning. A concrete defensive pattern:

def _care_solve(A, B, Q, R, N=None):
    from scipy.linalg import solve_continuous_are
    A_arr = np.asarray(A, dtype=float)
    B_arr = np.asarray(B, dtype=float)
    Q_arr = np.asarray(Q, dtype=float)
    R_arr = np.asarray(R, dtype=float)
    try:
        P = solve_continuous_are(A_arr, B_arr, Q_arr, R_arr, s=N)
    except np.linalg.LinAlgError as e:
        raise RuntimeError(f"CARE failed — check that (A, B) is stabilizable "
                           f"and (A, Q^{{1/2}}) is detectable. Underlying: {e}")
    # Defensive symmetry + positive-semidefiniteness check
    P = 0.5 * (P + P.T)  # symmetrize (CARE solution is mathematically symmetric)
    eig_min = np.linalg.eigvalsh(P).min()
    if eig_min < -1e-8 * np.linalg.norm(P):
        raise RuntimeError(f"CARE returned non-PSD P (min eig = {eig_min:.3e}). "
                           f"This usually means the Riccati solver failed silently.")
    return P

The symmetrization step is particularly important because floating-point accumulation in the Schur solver introduces asymmetry at the 1e-14 level that is harmless mathematically but propagates into the subsequent Cholesky solve for $K$.

3. Discretization at high sample rates — use zoh by default, not tustin.

The API proposal lists tustin as the default discretization method. Tustin (bilinear transform) is exact for sinusoidal inputs but introduces a frequency warp that becomes significant when $\omega_N / (\omega_s / 2) > 0.1$ (where $\omega_N$ is the natural frequency and $\omega_s / 2$ is the Nyquist frequency). In control-system practice, the sample rate is usually chosen to be 10-30x the dominant closed-loop bandwidth, which puts Tustin's warp error right in the regime where it starts to matter.

Zero-order hold (ZOH) is the correct default because it models the physical reality of a DAC driving the plant: the control input is held constant between samples. The formulas are:

$$ A_d = e^{A T_s}, \quad B_d = \left( \int_0^{T_s} e^{A \tau} d\tau \right) B $$

where $T_s$ is the sample period. scipy's cont2discrete(..., method='zoh') handles the matrix exponential correctly via Padé approximation. The behavior matches the physical system exactly under the held-input assumption, which is the correct assumption for any real digital control loop.

Tustin should still be available as an option for users who want it (e.g., for filter discretization where held-input doesn't make physical sense), but it should not be the default. I would make the default "zoh" and keep "tustin" as an explicit opt-in.

4. Stability margin computation — watch out for the multi-loop case.

The proposed gain_margin and phase_margin signatures take num and den arguments, which implies single-input single-output (SISO) systems. For a state-space system (MIMO in general), the classical gain/phase margin definitions do not apply directly — you need to compute the structured singular value (μ-analysis) or the disk margin (which is a newer SISO-like summary that extends cleanly to MIMO).

For the first version, I would explicitly scope the margin functions to SISO and either:

• Require the user to extract the SISO transfer function first (which is what the num, den signature already implies)
• Add a loop: int = 0 argument that lets the user specify which input-output channel to analyze
• Return None with a clear error message for MIMO systems

And add a TODO to implement disk margin or μ for v2 — it is ~100 lines on top of scipy's existing linear algebra and is what modern practice uses for robust control margins.

5. Gain scheduling as "LQR evaluated on a grid" is the right abstraction.

The issue mentions "gain scheduling" but the proposal does not flesh out the API. The standard approach is to:

1. Linearize the nonlinear plant at a set of operating points
2. Design an LQR gain $K_i$ for each linearization
3. Interpolate the gains $K(\theta) = \sum_i w_i(\theta) K_i$ as the scheduling variable $\theta$ varies

The tricky part is the interpolation: naive linear interpolation of the gain matrices can produce a scheduled closed-loop system that is unstable at the boundaries between grid cells even when the per-grid-cell designs are all stable. The fix is to interpolate the underlying Riccati solutions $P_i$ rather than the gains $K_i$, then derive the scheduled gain from the interpolated $P$:

$$ P(\theta) = \sum_i w_i(\theta) P_i, \quad K(\theta) = R^{-1} B(\theta)^T P(\theta) $$

This is a paper-level result from Stilwell & Rugh 2000 that is not widely implemented, and it would be a genuine differentiator for gds-analysis over the usual "just interpolate gains" approach. Worth at least documenting in the API even if the first implementation is the naive version.

6. One small API consistency note.

The proposal uses list[list[float]] as the universal matrix type to avoid LinearizedSystem coupling, which is the right choice for decoupling. One thing to watch: when returning complex eigenvalues (e.g., from eigenvalues(A)), list[complex] is the natural choice, but JSON-serialization-friendly code often wants list[tuple[float, float]] (real, imag). Consider offering both via a as_complex_tuple: bool = False flag, or document clearly that the return is Python complex.

Happy to dig deeper on any of these once the PR is up — the CARE solver's ill-conditioning and the gain-scheduling Stilwell-Rugh trick are the two points I would want to get into the initial implementation rather than left for v2.