From c04715e5b2ca32894eb9d6d99a1a797f7b474278 Mon Sep 17 00:00:00 2001 From: rohan Date: Tue, 7 Apr 2026 18:44:31 +0530 Subject: [PATCH 1/3] feat: add classical control theory analysis stack (#198) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add six new modules across four packages implementing the control theory capabilities identified in the MATLAB tech talk gap analysis: gds-domains/symbolic: - transfer.py: transfer functions (ss_to_tf), poles/zeros, controllability/observability, Gang of Six sensitivity (#199) - delay.py: Padé approximation for time delay modeling (#200) gds-analysis: - linear.py: eigenvalue stability, frequency response, gain/phase margins, discretization (Tustin/ZOH/Euler), LQR/DLQR, Kalman filter, gain scheduling (#201) - response.py: step/impulse response computation and time-domain metrics (rise time, settling time, overshoot, SSE) (#202) gds-proof/analysis: - lyapunov.py: Lyapunov candidate verification, quadratic Lyapunov (A'P+PA), Lyapunov equation solver, passivity certificates (#203) gds-viz: - frequency.py: Bode, Nyquist, Nichols, root locus plots (#204) - response.py: annotated step/impulse response plots Zero new third-party dependencies — all capabilities fit within existing optional extras ([symbolic], [continuous], [phase]/[control]). 78 new tests, all existing tests unaffected. --- improvement-plans/review_synthesis.md | 166 ------ .../gds-analysis/gds_analysis/__init__.py | 44 +- packages/gds-analysis/gds_analysis/linear.py | 557 ++++++++++++++++++ .../gds-analysis/gds_analysis/response.py | 328 +++++++++++ packages/gds-analysis/tests/test_linear.py | 236 ++++++++ packages/gds-analysis/tests/test_response.py | 152 +++++ .../gds_domains/symbolic/__init__.py | 26 + .../gds-domains/gds_domains/symbolic/delay.py | 120 ++++ .../gds_domains/symbolic/transfer.py | 437 ++++++++++++++ packages/gds-domains/tests/test_delay.py | 99 ++++ packages/gds-domains/tests/test_transfer.py | 302 ++++++++++ .../gds-proof/gds_proof/analysis/lyapunov.py | 527 +++++++++++++++++ packages/gds-proof/tests/test_lyapunov.py | 175 ++++++ packages/gds-viz/gds_viz/__init__.py | 22 +- packages/gds-viz/gds_viz/frequency.py | 455 ++++++++++++++ packages/gds-viz/gds_viz/response.py | 272 +++++++++ packages/gds-viz/pyproject.toml | 2 + 17 files changed, 3750 insertions(+), 170 deletions(-) delete mode 100644 improvement-plans/review_synthesis.md create mode 100644 packages/gds-analysis/gds_analysis/linear.py create mode 100644 packages/gds-analysis/gds_analysis/response.py create mode 100644 packages/gds-analysis/tests/test_linear.py create mode 100644 packages/gds-analysis/tests/test_response.py create mode 100644 packages/gds-domains/gds_domains/symbolic/delay.py create mode 100644 packages/gds-domains/gds_domains/symbolic/transfer.py create mode 100644 packages/gds-domains/tests/test_delay.py create mode 100644 packages/gds-domains/tests/test_transfer.py create mode 100644 packages/gds-proof/gds_proof/analysis/lyapunov.py create mode 100644 packages/gds-proof/tests/test_lyapunov.py create mode 100644 packages/gds-viz/gds_viz/frequency.py create mode 100644 packages/gds-viz/gds_viz/response.py diff --git a/improvement-plans/review_synthesis.md b/improvement-plans/review_synthesis.md deleted file mode 100644 index 225cec9..0000000 --- a/improvement-plans/review_synthesis.md +++ /dev/null @@ -1,166 +0,0 @@ -# Synthesis of GDS-Core Reviews: Zargham vs. Jamsheed - -**Date:** 2026-04-03 -**Sources:** -- `gds_core_briefing.md`, `gds_core_improvement_roadmap.md`, `gds_core_risk_register_and_doctrines.md` — Engineering review (Zargham + Claude understudy) -- `math-review.md` — Mathematical review (Jamsheed) - ---- - -## Executive Summary - -Both reviews validate the theoretical foundation but identify critical gaps in notation alignment with the paper, `ControlAction` semantics, and temporal agnosticism documentation. Jamsheed focuses on mathematical formalism fidelity to Zargham & Shorish (2022); Zargham provides a 14-item engineering roadmap with risks and doctrines. Their recommendations are ~85% aligned, with minor divergences on tier ordering. - ---- - -## Where the Reviews Agree - -### 1. ControlAction = Output Map (Not Admissibility Constraint) - -| Review | Finding | -|--------|---------| -| **Zargham T0-3** | `ControlAction` is unused by all six DSLs; should be `y = C(x,d)` for controller-plant duality | -| **Jamsheed** | "erroneously referred to in the codebase as 'admissibility constraint'" — should be output observable for composition | - -**Consensus:** Same signal `y` is: -- **Inside perspective:** plant's output map -- **Outside perspective:** control action on the next system - -### 2. Notation Must Harmonize with Paper - -| Issue | Codebase | Paper/Bosch | Fix | -|-------|----------|-------------|-----| -| External factors | `u` | `z` (or `varepsilon`) | Rename to `z` or distinct from action `u` | -| "Input Space U" | Used for external factors | `U_x` is admissible input space | Clarify terminology | -| Policy mapping | `g(x,u)` | `g(x)` or `g(x,z)` | Align with paper's intent | - -### 3. Core is Time-Agnostic (Structural Recurrence) - -| Review | Framing | -|--------|---------| -| **Zargham T0-4** | "Core is time-agnostic" — `is_temporal=True` is structural boundary marker only | -| **Jamsheed** | "Structural recurrence primitive regardless of `is_temporal`" — uses semigroup `(N, +)` for sequencing | - -**Consensus:** The three-layer temporal stack: -- Layer 0 (core): Structural recurrence only — no time model -- Layer 1 (DSL): `ExecutionContract` declares time model -- Layer 2 (sim): Solver/runner implements time model - -### 4. Representation-Agnostic Architecture - -Both confirm the core should not commit to numerical, relational, or any substrate. Arithmetic enters at DSL layer. - ---- - -## Where the Reviews Diverge - -### Tier Ordering Disagreements - -| Item | Zargham's Tier | Jamsheed's View | -|------|---------------|-----------------| -| **Continuous-time formalization** | Tier 2 (medium priority) | Move to Tier 1 — "ambitious objective... many limit conceptions possible" | -| **Stochastic extensions** | Tier 3 (research frontier) | "Most crucial... could be Tier in own right or Tier 2" | -| **Discrete-time framework** | Moore as default in T1-1 | Elevate Mealy machines — "Mealy includes Moore as subclass" | - -### Scope Differences - -| Aspect | Zargham | Jamsheed | -|--------|---------|----------| -| **Primary output** | 14-item roadmap with risks/doctrines | Mathematical corrections + architectural feedback | -| **ControlAction taxonomy** | Deep dive on controller-plant duality, perspective inversion, observability connection | Suggests `Output`/`OutputObservable` class with typing for plant vs control action | -| **Time agnosticism** | Three-layer stack with formal invariant | Emphasizes "before/after" semantic structure; notes OGS proves degenerate case works | - ---- - -## The Criticisms: A Taxonomy - -### Tier 0 (Must Fix Before Anything Else) - -| ID | Issue | Severity | Reviews Citing | -|----|-------|----------|---------------| -| T0-1 | Checks lack formal property statements | **Critical** | Zargham (C1 violation) | -| T0-2 | Tests lack requirement traceability | **Critical** | Zargham (C1, C3) | -| T0-3 | `ControlAction` role unresolved | **Critical** | **Both** — Zargham's T0-3, Jamsheed's General Comments | -| T0-4 | Temporal agnosticism not formally stated | **Critical** | **Both** — Zargham's T0-4, Jamsheed's time agnosticism section | - -### Documentation Issues - -| Issue | Jamsheed's Words | Zargham's Mitigation | -|-------|------------------|---------------------| -| `ControlAction` misnamed | "erroneously referred to... as admissibility constraint" | T0-3 deliverable 2: dedicated duality documentation page | -| Input space conflation | "erroneously referred to as the 'Input Space U'" | OC-6: distinct categories for `U_c`, `W`, `D`, `Y` | -| Discrete-time bias | "documentation quietly undermines [time agnosticism]" | T0-4 deliverable 3: audit + correction of trajectory notation | - -### Ontological Commitments - -Jamsheed's review confirms Zargham's OC-1 through OC-10 but adds nuance: - -> "Generalized refers to general data structures (representation generality, architecture principle 1), while dynamical refers to both the trajectory sequencing and its computational implementability (architecture principles 2 and 3)." - -This directly validates Zargham's: -- OC-3 (Space as compatibility token — "generalized") -- OC-7 (Structural recurrence — "dynamical") -- OC-10 (Simulation is instrumentation — "computational implementability") - ---- - -## Actionable Steps: Unified Priority - -### Immediate (Tier 0 — No New Code) - -**1. T0-3: Resolve ControlAction** — Both reviews agree this is foundational -- Formalize `y = C(x,d)` as output map -- Document controller-plant duality with thermostat example -- Add `(Y, C)` to canonical projection -- Create SC-check preventing ControlAction in `g` pathway - -**2. Notation Harmonization** — Jamsheed's math review + Zargham's OC-6 -- External factors: change `u` to `z` (or distinct symbol) -- Rename "Input Space U" to "External Factor Space Z" -- Policy: clarify `g(x)` vs `g(x,z)` vs `g(x,u_c)` with controlled inputs - -**3. T0-4: Formalize Temporal Agnosticism** -- State invariant: `is_temporal=True` = structural recurrence only -- Three-layer stack diagram -- Audit docs: replace "next step" with "temporal boundary" - -### Near-Term (Tier 1) - -**4. T1-1: ExecutionContract** — Elevate per Jamsheed's feedback on Mealy machines -- Add `update_ordering: Literal["Moore", "Mealy"]` field -- Default remains Moore, but Mealy is explicitly supported - -**5. T1-3: Disturbance Inputs** — Resolves Jamsheed's `z` notation concern -- Formal partition: `U_c` (controlled) vs `W` (disturbance) -- Tagged mixin: `{"role": "disturbance"}` -- DST-001 check: disturbance bypasses policy layer - -### Medium-Term (Tier 2 — Per Jamsheed's Recommendations) - -**6. T2-4: Continuous-Time** — Consider moving to Tier 1 per Jamsheed -- `SolverInterface` contract separates spec from simulation -- Scope: real-valued subspaces only (explicit constraint) - -**7. Stochastic Extensions** — Jamsheed's "most crucial" item -- Consider Tier 2 placement (not Tier 3) per Jamsheed's recommendation -- Requires T1-1 (`ExecutionContract`) and T2-2 (behavioral verification) stable first - ---- - -## Key Architectural Decisions Pending - -1. **Should continuous-time formalization move to Tier 1?** Jamsheed argues yes for earlier operationalization; Zargham placed it Tier 2 for scope discipline. - -2. **Should stochastic extensions be elevated to Tier 2 (or own Tier)?** Jamsheed sees this as crucial; roadmap currently has it Tier 3. - -3. **Mealy vs Moore as discrete-time framework:** Jamsheed suggests elevating Mealy as the general case (Moore as subclass). Currently roadmap has Moore as default with Mealy as option. - -4. **External factor notation:** `z` (Bosch presentation), `varepsilon` (past usage), or keep `u` with clear subscripting? Jamsheed flags conflict with action `u`. - ---- - -## Changelog - -| Version | Date | Change | -|---------|------|--------| -| 0.1 | 2026-04-03 | Initial synthesis from both reviews | diff --git a/packages/gds-analysis/gds_analysis/__init__.py b/packages/gds-analysis/gds_analysis/__init__.py index 2e561e4..6ef07eb 100644 --- a/packages/gds-analysis/gds_analysis/__init__.py +++ b/packages/gds-analysis/gds_analysis/__init__.py @@ -1,8 +1,8 @@ """Dynamical analysis for GDS specifications. Bridges gds-framework structural annotations to gds-sim runtime, -enabling constraint enforcement, metric computation, and reachability -analysis on concrete trajectories. +enabling constraint enforcement, metric computation, reachability +analysis, linear systems analysis, and response metrics. """ __version__ = "0.1.2" @@ -22,11 +22,13 @@ reachable_graph, reachable_set, ) +from gds_analysis.response import StepMetrics, step_response_metrics __all__ = [ "BackwardReachableSet", "Isochrone", "ReachabilityResult", + "StepMetrics", "backward_reachable_set", "configuration_space", "extract_isochrones", @@ -34,5 +36,43 @@ "reachable_graph", "reachable_set", "spec_to_model", + "step_response_metrics", "trajectory_distances", ] + + +# Lazy imports for optional modules (require [continuous] extra) +_LINEAR_EXPORTS = { + "eigenvalues", + "is_stable", + "is_marginally_stable", + "frequency_response", + "gain_margin", + "phase_margin", + "discretize", + "lqr", + "dlqr", + "kalman", + "gain_schedule", +} + +_RESPONSE_SCIPY_EXPORTS = { + "step_response", + "impulse_response", + "metrics_from_ode_results", +} + + +def __getattr__(name: str) -> object: + """Lazy import for optional scipy-dependent functions.""" + if name in _LINEAR_EXPORTS: + from gds_analysis import linear + + return getattr(linear, name) + + if name in _RESPONSE_SCIPY_EXPORTS: + from gds_analysis import response + + return getattr(response, name) + + raise AttributeError(f"module 'gds_analysis' has no attribute {name!r}") diff --git a/packages/gds-analysis/gds_analysis/linear.py b/packages/gds-analysis/gds_analysis/linear.py new file mode 100644 index 0000000..95d7468 --- /dev/null +++ b/packages/gds-analysis/gds_analysis/linear.py @@ -0,0 +1,557 @@ +"""Numerical linear systems analysis, discretization, and controller synthesis. + +Provides eigenvalue stability checks, frequency response computation, +stability margins, continuous-to-discrete conversion, LQR/Kalman synthesis, +and gain scheduling. + +All functions accept ``list[list[float]]`` matrices (matching +``LinearizedSystem`` field types from gds-domains/symbolic) and return +plain Python types. Numpy/scipy are internal implementation details. + +Requires ``gds-analysis[continuous]`` (scipy + numpy). +""" + +from __future__ import annotations + +from typing import Any + +_IMPORT_MSG = ( + "gds_analysis.linear requires scipy and numpy. " + "Install with: uv add gds-analysis[continuous]" +) + + +def _require_deps() -> None: + """Raise ImportError if scipy/numpy are absent.""" + try: + import numpy # noqa: F401 + import scipy # noqa: F401 + except ImportError as exc: + raise ImportError(_IMPORT_MSG) from exc + + +# --------------------------------------------------------------------------- +# Stability analysis +# --------------------------------------------------------------------------- + + +def eigenvalues(A: list[list[float]]) -> list[complex]: + """Eigenvalues of state matrix A. + + Parameters + ---------- + A : list[list[float]] + Square state matrix. + + Returns + ------- + list[complex] + Eigenvalues (may have zero imaginary part for real eigenvalues). + """ + _require_deps() + import numpy as np + + if not A: + return [] + return [complex(e) for e in np.linalg.eigvals(np.array(A, dtype=float))] + + +def is_stable(A: list[list[float]], *, continuous: bool = True) -> bool: + """Check asymptotic stability of the state matrix. + + Parameters + ---------- + A : list[list[float]] + Square state matrix. + continuous : bool + If True (default), checks continuous-time stability: all Re(λ) < 0. + If False, checks discrete-time stability: all |λ| < 1. + + Returns + ------- + bool + True if the system is asymptotically stable. + """ + eigs = eigenvalues(A) + if not eigs: + return True # vacuously stable + + if continuous: + return all(e.real < 0 for e in eigs) + return all(abs(e) < 1.0 for e in eigs) + + +def is_marginally_stable(A: list[list[float]], *, continuous: bool = True) -> bool: + """Check marginal stability. + + A system is marginally stable if no eigenvalue is in the unstable + region and at least one is on the stability boundary. + + Parameters + ---------- + A : list[list[float]] + Square state matrix. + continuous : bool + If True, boundary is Re(λ) = 0. If False, boundary is |λ| = 1. + """ + eigs = eigenvalues(A) + if not eigs: + return False + + tol = 1e-10 + if continuous: + if any(e.real > tol for e in eigs): + return False + return any(abs(e.real) <= tol for e in eigs) + else: + if any(abs(e) > 1.0 + tol for e in eigs): + return False + return any(abs(abs(e) - 1.0) <= tol for e in eigs) + + +# --------------------------------------------------------------------------- +# 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, + *, + n_points: int = 500, + omega_range: tuple[float, float] = (1e-2, 1e2), +) -> tuple[list[float], list[float], list[float]]: + """Compute frequency response H(jω) numerically. + + Parameters + ---------- + A, B, C, D : list[list[float]] + State-space matrices (SISO or first input/output for MIMO). + omega : list[float] | None + Frequency points (rad/s). If None, auto-generates from omega_range. + n_points : int + Number of frequency points if omega is None. + omega_range : tuple[float, float] + Log-spaced frequency range if omega is None. + + Returns + ------- + (omega, magnitude_db, phase_deg) + All as plain ``list[float]``. + """ + _require_deps() + import numpy as np + from scipy import signal + + sys = signal.StateSpace( + np.array(A, dtype=float), + np.array(B, dtype=float), + np.array(C, dtype=float), + np.array(D, dtype=float), + ) + + if omega is None: + w = np.logspace(np.log10(omega_range[0]), np.log10(omega_range[1]), n_points) + else: + w = np.array(omega, dtype=float) + + w_out, H = signal.freqresp(sys, w) + + # For MIMO, take [0,0] element + if H.ndim > 1: + H = H[0, 0, :] + + mag_db = (20 * np.log10(np.abs(H))).tolist() + phase_deg = (np.degrees(np.unwrap(np.angle(H)))).tolist() + + return w_out.tolist(), mag_db, phase_deg + + +def gain_margin(num: list[float], den: list[float]) -> tuple[float, float]: + """Compute gain margin. + + Parameters + ---------- + num : list[float] + Transfer function numerator coefficients (descending powers). + den : list[float] + Transfer function denominator coefficients (descending powers). + + Returns + ------- + (gm_db, gm_freq) + Gain margin in dB and the frequency (rad/s) where phase = -180°. + Returns (float('inf'), float('nan')) if phase never crosses -180°. + """ + _require_deps() + import numpy as np + from scipy import signal + + sys = signal.TransferFunction( + np.array(num, dtype=float), + np.array(den, dtype=float), + ) + w = np.logspace(-4, 4, 10000) + w_out, H = signal.freqresp(sys, w) + + phase = np.unwrap(np.angle(H)) + + # Find -180° crossing + target = -np.pi + crossings = [] + for i in range(1, len(phase)): + if (phase[i - 1] - target) * (phase[i] - target) < 0: + # Interpolate + frac = (target - phase[i - 1]) / (phase[i] - phase[i - 1]) + w_cross = w_out[i - 1] + frac * (w_out[i] - w_out[i - 1]) + mag_at_cross = np.abs(H[i - 1]) + frac * (np.abs(H[i]) - np.abs(H[i - 1])) + crossings.append((w_cross, mag_at_cross)) + + if not crossings: + return float("inf"), float("nan") + + # Take the first crossing + w_cross, mag_at_cross = crossings[0] + gm_db = float(-20 * np.log10(mag_at_cross)) if mag_at_cross > 0 else float("inf") + + return gm_db, float(w_cross) + + +def phase_margin(num: list[float], den: list[float]) -> tuple[float, float]: + """Compute phase margin. + + Parameters + ---------- + num : list[float] + Transfer function numerator coefficients (descending powers). + den : list[float] + Transfer function denominator coefficients (descending powers). + + Returns + ------- + (pm_deg, pm_freq) + Phase margin in degrees and the gain crossover frequency (rad/s). + Returns (float('inf'), float('nan')) if gain never crosses 0 dB. + """ + _require_deps() + import numpy as np + from scipy import signal + + sys = signal.TransferFunction( + np.array(num, dtype=float), + np.array(den, dtype=float), + ) + w = np.logspace(-4, 4, 10000) + w_out, H = signal.freqresp(sys, w) + + mag = np.abs(H) + phase = np.unwrap(np.angle(H)) + + # Find 0 dB crossing (|H| = 1) + crossings = [] + for i in range(1, len(mag)): + if (mag[i - 1] - 1.0) * (mag[i] - 1.0) < 0: + frac = (1.0 - mag[i - 1]) / (mag[i] - mag[i - 1]) + w_cross = w_out[i - 1] + frac * (w_out[i] - w_out[i - 1]) + phase_at_cross = phase[i - 1] + frac * (phase[i] - phase[i - 1]) + crossings.append((w_cross, phase_at_cross)) + + if not crossings: + return float("inf"), float("nan") + + # Take the first crossing + w_cross, phase_at_cross = crossings[0] + pm_deg = float(np.degrees(phase_at_cross) + 180.0) + + return pm_deg, float(w_cross) + + +# --------------------------------------------------------------------------- +# 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). + + Parameters + ---------- + A, B, C, D : list[list[float]] + Continuous-time state-space matrices. + dt : float + Sampling period in seconds. + method : str + Discretization method: ``"tustin"`` (bilinear), ``"zoh"`` + (zero-order hold), ``"euler"`` (forward Euler). + + Returns + ------- + (Ad, Bd, Cd, Dd) as ``list[list[float]]``. + """ + _require_deps() + import numpy as np + from scipy import signal + + method_map = { + "tustin": "bilinear", + "bilinear": "bilinear", + "zoh": "zoh", + "euler": "euler", + "forward_euler": "euler", + "backward_euler": "backward_diff", + "backward_diff": "backward_diff", + } + + scipy_method = method_map.get(method.lower()) + if scipy_method is None: + raise ValueError( + f"Unknown discretization method '{method}'. " + f"Valid methods: {sorted(method_map.keys())}" + ) + + sys_c = signal.StateSpace( + np.array(A, dtype=float), + np.array(B, dtype=float), + np.array(C, dtype=float), + np.array(D, dtype=float), + ) + + sys_d = sys_c.to_discrete(dt, method=scipy_method) + + return ( + sys_d.A.tolist(), + sys_d.B.tolist(), + sys_d.C.tolist(), + sys_d.D.tolist(), + ) + + +# --------------------------------------------------------------------------- +# 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]]: + """Continuous-time Linear Quadratic Regulator. + + Minimizes J = ∫(x'Qx + u'Ru + 2x'Nu) dt by solving the continuous + algebraic Riccati equation. + + Parameters + ---------- + A : list[list[float]] + State matrix (n x n). + B : list[list[float]] + Input matrix (n x m). + Q : list[list[float]] + State cost matrix (n x n, positive semi-definite). + R : list[list[float]] + Input cost matrix (m x m, positive definite). + N : list[list[float]] | None + Cross-term matrix (n x m). Default is zero. + + Returns + ------- + (K, P, E) + K — optimal gain matrix (u = -Kx), shape (m x n). + P — solution to the Riccati equation, shape (n x n). + E — closed-loop eigenvalues. + + Raises + ------ + ValueError + If the system is not stabilizable. + """ + _require_deps() + import numpy as np + from scipy import linalg + + a = np.array(A, dtype=float) + b = np.array(B, dtype=float) + q = np.array(Q, dtype=float) + r = np.array(R, dtype=float) + + if N is not None: + n_mat = np.array(N, dtype=float) + # Transform: A_bar = A - B R^{-1} N', Q_bar = Q - N R^{-1} N' + r_inv = np.linalg.inv(r) + a_bar = a - b @ r_inv @ n_mat.T + q_bar = q - n_mat @ r_inv @ n_mat.T + else: + a_bar = a + q_bar = q + n_mat = np.zeros((a.shape[0], b.shape[1])) + + try: + P = linalg.solve_continuous_are(a_bar, b, q_bar, r) + except np.linalg.LinAlgError as exc: + raise ValueError( + "Failed to solve Riccati equation. " + "Check that (A, B) is stabilizable and Q >= 0, R > 0." + ) from exc + + r_inv = np.linalg.inv(r) + K = r_inv @ (b.T @ P + n_mat.T) + + # Closed-loop eigenvalues + A_cl = a - b @ K + E = np.linalg.eigvals(A_cl) + + return K.tolist(), P.tolist(), [complex(e) for e in E] + + +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 Linear Quadratic Regulator. + + Minimizes J = Σ(x'Qx + u'Ru) by solving the discrete algebraic + Riccati equation. + + Parameters + ---------- + Ad : list[list[float]] + Discrete state matrix (n x n). + Bd : list[list[float]] + Discrete input matrix (n x m). + Q : list[list[float]] + State cost (n x n, positive semi-definite). + R : list[list[float]] + Input cost (m x m, positive definite). + + Returns + ------- + (K, P, E) + K — optimal gain, P — Riccati solution, E — closed-loop eigenvalues. + """ + _require_deps() + import numpy as np + from scipy import linalg + + a = np.array(Ad, dtype=float) + b = np.array(Bd, dtype=float) + q = np.array(Q, dtype=float) + r = np.array(R, dtype=float) + + try: + P = linalg.solve_discrete_are(a, b, q, r) + except np.linalg.LinAlgError as exc: + raise ValueError( + "Failed to solve discrete Riccati equation. " + "Check that (Ad, Bd) is stabilizable." + ) from exc + + K = np.linalg.inv(r + b.T @ P @ b) @ (b.T @ P @ a) + + A_cl = a - b @ K + E = np.linalg.eigvals(A_cl) + + return K.tolist(), P.tolist(), [complex(e) for e in E] + + +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 the continuous algebraic Riccati equation for + the observer problem. + + Parameters + ---------- + A : list[list[float]] + State matrix (n x n). + C : list[list[float]] + Output matrix (p x n). + Q_process : list[list[float]] + Process noise covariance (n x n). + R_measurement : list[list[float]] + Measurement noise covariance (p x p). + + Returns + ------- + (L, P) + L — observer gain matrix (n x p), such that x̂' = Ax̂ + Bu + L(y - Cx̂). + P — error covariance matrix (n x n). + + Raises + ------ + ValueError + If (A, C) is not detectable. + """ + _require_deps() + import numpy as np + from scipy import linalg + + a = np.array(A, dtype=float) + c = np.array(C, dtype=float) + q = np.array(Q_process, dtype=float) + r = np.array(R_measurement, dtype=float) + + # Dual of LQR: solve ARE for A', C', Q, R + try: + P = linalg.solve_continuous_are(a.T, c.T, q, r) + except np.linalg.LinAlgError as exc: + raise ValueError( + "Failed to solve observer Riccati equation. " + "Check that (A, C) is detectable." + ) from exc + + L = P @ c.T @ np.linalg.inv(r) + + return L.tolist(), P.tolist() + + +def gain_schedule( + linearize_fn: Any, + operating_points: list[dict[str, float]], + Q: list[list[float]], + R: list[list[float]], +) -> list[tuple[dict[str, float], list[list[float]], list[complex]]]: + """Compute LQR gains at multiple operating points. + + Parameters + ---------- + linearize_fn : callable + Takes an operating point dict, returns + ``(A, B, C, D)`` as a tuple of ``list[list[float]]``. + operating_points : list[dict] + Operating point parameter dicts. + Q : list[list[float]] + State cost matrix (same for all points). + R : list[list[float]] + Input cost matrix (same for all points). + + Returns + ------- + list of (operating_point, K, closed_loop_eigenvalues) tuples. + """ + results = [] + for point in operating_points: + A_op, B_op, _C_op, _D_op = linearize_fn(point) + K, _P, E = lqr(A_op, B_op, Q, R) + results.append((point, K, E)) + return results diff --git a/packages/gds-analysis/gds_analysis/response.py b/packages/gds-analysis/gds_analysis/response.py new file mode 100644 index 0000000..517a93d --- /dev/null +++ b/packages/gds-analysis/gds_analysis/response.py @@ -0,0 +1,328 @@ +"""Step and impulse response computation with time-domain metrics. + +Provides step/impulse response generation from state-space models and +extraction of standard performance metrics (rise time, settling time, +overshoot, steady-state error). + +Requires ``gds-analysis[continuous]`` for response generation. +Metrics computation (``step_response_metrics``) works on any time-series +data and has no optional dependencies. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any + + +@dataclass(frozen=True) +class StepMetrics: + """Standard time-domain performance metrics for a step response. + + All time values are in the same units as the input time array. + """ + + rise_time: float + """Time from ``rise_lower`` to ``rise_upper`` fraction of final value.""" + + settling_time: float + """Time to permanently enter the ±settling_band around final value.""" + + overshoot_pct: float + """Peak overshoot as a percentage of the final value. Zero if none.""" + + peak_time: float + """Time at which the first peak occurs.""" + + peak_value: float + """Response value at the first peak.""" + + steady_state_value: float + """Estimated final value (mean of last 5% of samples).""" + + steady_state_error: float + """Absolute difference between setpoint and steady-state value.""" + + +def step_response_metrics( + times: list[float], + values: list[float], + setpoint: float = 1.0, + *, + settling_band: float = 0.02, + rise_lower: float = 0.1, + rise_upper: float = 0.9, +) -> StepMetrics: + """Extract standard performance metrics from a step response. + + Works on any time-series data — can be fed from: + + - ``step_response()`` output (this module) + - ``ODEResults.state_array()`` (gds-continuous) + - gds-sim ``Results`` column extraction + + Parameters + ---------- + times : list[float] + Monotonically increasing time values. + values : list[float] + Response values at each time point. + setpoint : float + Desired final value (default 1.0 for unit step). + settling_band : float + Fractional band for settling time (default ±2%). + rise_lower : float + Lower fraction for rise time (default 10%). + rise_upper : float + Upper fraction for rise time (default 90%). + + Returns + ------- + StepMetrics + Computed performance metrics. + """ + if len(times) < 2: + raise ValueError("Need at least 2 data points") + if len(times) != len(values): + raise ValueError("times and values must have the same length") + + n = len(values) + + # Steady-state value: mean of last 5% of samples + tail_count = max(1, n // 20) + ss_value = sum(values[-tail_count:]) / tail_count + ss_error = abs(setpoint - ss_value) + + # Initial value + v0 = values[0] + + # Range from initial to steady-state + v_range = ss_value - v0 + if abs(v_range) < 1e-15: + # No change — return degenerate metrics + return StepMetrics( + rise_time=0.0, + settling_time=0.0, + overshoot_pct=0.0, + peak_time=times[0], + peak_value=values[0], + steady_state_value=ss_value, + steady_state_error=ss_error, + ) + + # Rise time: time from rise_lower to rise_upper fraction of final value + low_threshold = v0 + rise_lower * v_range + high_threshold = v0 + rise_upper * v_range + + t_low = _interpolate_crossing(times, values, low_threshold) + t_high = _interpolate_crossing(times, values, high_threshold) + rise_time = t_high - t_low if t_high > t_low else 0.0 + + # Peak value and time + if v_range > 0: + peak_idx = max(range(n), key=lambda i: values[i]) + else: + peak_idx = min(range(n), key=lambda i: values[i]) + peak_value = values[peak_idx] + peak_time = times[peak_idx] + + # Overshoot percentage + overshoot = (peak_value - ss_value) / abs(v_range) * 100.0 + overshoot_pct = max(0.0, overshoot) if v_range > 0 else max(0.0, -overshoot) + + # Settling time: backward scan + band = settling_band * abs(ss_value) if abs(ss_value) > 1e-15 else settling_band + settling_time = times[0] + for i in range(n - 1, -1, -1): + if abs(values[i] - ss_value) > band: + settling_time = times[min(i + 1, n - 1)] + break + + return StepMetrics( + rise_time=rise_time, + settling_time=settling_time, + overshoot_pct=overshoot_pct, + peak_time=peak_time, + peak_value=peak_value, + steady_state_value=ss_value, + steady_state_error=ss_error, + ) + + +def _interpolate_crossing( + times: list[float], + values: list[float], + threshold: float, +) -> float: + """Find the first time at which values crosses threshold, with interpolation.""" + for i in range(1, len(values)): + v_prev, v_curr = values[i - 1], values[i] + if (v_prev <= threshold <= v_curr) or (v_prev >= threshold >= v_curr): + # Linear interpolation + dv = v_curr - v_prev + if abs(dv) < 1e-15: + return times[i] + frac = (threshold - v_prev) / dv + return times[i - 1] + frac * (times[i] - times[i - 1]) + # Threshold never crossed — return last time + return times[-1] + + +def step_response( + A: list[list[float]], + B: list[list[float]], + C: list[list[float]], + D: list[list[float]], + t_span: tuple[float, float] = (0.0, 10.0), + n_points: int = 1000, + *, + input_index: int = 0, +) -> tuple[list[float], list[list[float]]]: + """Compute unit step response from a state-space model. + + Applies a unit step to the specified input and simulates. + + Parameters + ---------- + A, B, C, D : list[list[float]] + State-space matrices. + t_span : tuple[float, float] + Simulation time interval. + n_points : int + Number of evaluation points. + input_index : int + Which input channel receives the step (default 0). + + Returns + ------- + (times, outputs) + ``times`` is a list of float. ``outputs[i]`` is the response + of output i as a list of float. + + Requires + -------- + ``gds-analysis[continuous]`` (scipy + numpy). + """ + try: + import numpy as np + from scipy import signal + except ImportError as exc: + raise ImportError( + "step_response requires scipy and numpy. " + "Install with: uv add gds-analysis[continuous]" + ) from exc + + sys = signal.StateSpace( + np.array(A, dtype=float), + np.array(B, dtype=float), + np.array(C, dtype=float), + np.array(D, dtype=float), + ) + t_eval = np.linspace(t_span[0], t_span[1], n_points) + + # Use lsim with step input + n_inputs = np.array(B).shape[1] if len(A) > 0 else np.array(D).shape[1] + u_input = np.zeros((n_points, n_inputs)) + u_input[:, input_index] = 1.0 + + tout, yout, _ = signal.lsim(sys, U=u_input, T=t_eval) + + times = tout.tolist() + if yout.ndim == 1: + outputs = [yout.tolist()] + else: + outputs = [yout[:, i].tolist() for i in range(yout.shape[1])] + + return times, outputs + + +def impulse_response( + A: list[list[float]], + B: list[list[float]], + C: list[list[float]], + D: list[list[float]], + t_span: tuple[float, float] = (0.0, 10.0), + n_points: int = 1000, + *, + input_index: int = 0, +) -> tuple[list[float], list[list[float]]]: + """Compute impulse response from a state-space model. + + Parameters + ---------- + A, B, C, D : list[list[float]] + State-space matrices. + t_span : tuple[float, float] + Simulation time interval. + n_points : int + Number of evaluation points. + input_index : int + Which input channel receives the impulse (default 0). + + Returns + ------- + (times, outputs) + ``times`` is a list of float. ``outputs[i]`` is the response + of output i as a list of float. + + Requires + -------- + ``gds-analysis[continuous]`` (scipy + numpy). + """ + try: + import numpy as np + from scipy import signal + except ImportError as exc: + raise ImportError( + "impulse_response requires scipy and numpy. " + "Install with: uv add gds-analysis[continuous]" + ) from exc + + a = np.array(A, dtype=float) + b = np.array(B, dtype=float) + c = np.array(C, dtype=float) + d = np.array(D, dtype=float) + + # Extract SISO for the specified input + b_col = b[:, input_index : input_index + 1] if len(A) > 0 else np.array([[0.0]]) + d_col = d[:, input_index : input_index + 1] + + sys = signal.StateSpace(a, b_col, c, d_col) + t_eval = np.linspace(t_span[0], t_span[1], n_points) + tout, yout = signal.impulse(sys, T=t_eval) + + times = tout.tolist() + if yout.ndim == 1: + outputs = [yout.tolist()] + else: + outputs = [yout[:, i].tolist() for i in range(yout.shape[1])] + + return times, outputs + + +def metrics_from_ode_results( + results: Any, + state_name: str, + setpoint: float = 1.0, + **kwargs: Any, +) -> StepMetrics: + """Extract step metrics from an ``ODEResults`` object. + + Convenience wrapper: calls ``results.times`` and + ``results.state_array(state_name)``, then delegates to + ``step_response_metrics()``. + + Parameters + ---------- + results + An ``ODEResults`` object from gds-continuous. + state_name : str + Name of the state variable to analyze. + setpoint : float + Desired final value. + **kwargs + Passed through to ``step_response_metrics()``. + """ + times = results.times + values = results.state_array(state_name) + return step_response_metrics(times, values, setpoint, **kwargs) diff --git a/packages/gds-analysis/tests/test_linear.py b/packages/gds-analysis/tests/test_linear.py new file mode 100644 index 0000000..084ec34 --- /dev/null +++ b/packages/gds-analysis/tests/test_linear.py @@ -0,0 +1,236 @@ +"""Tests for numerical linear systems analysis.""" + +import pytest + +from gds_analysis.linear import ( + discretize, + dlqr, + eigenvalues, + frequency_response, + gain_margin, + is_marginally_stable, + is_stable, + kalman, + lqr, + phase_margin, +) + + +# --------------------------------------------------------------------------- +# Stability analysis +# --------------------------------------------------------------------------- + + +class TestEigenvalues: + def test_stable_diagonal(self) -> None: + eigs = eigenvalues([[-1.0, 0.0], [0.0, -2.0]]) + reals = sorted([e.real for e in eigs]) + assert reals[0] == pytest.approx(-2.0, abs=1e-10) + assert reals[1] == pytest.approx(-1.0, abs=1e-10) + + def test_unstable_system(self) -> None: + eigs = eigenvalues([[1.0, 0.0], [0.0, -1.0]]) + reals = sorted([e.real for e in eigs]) + assert reals[0] == pytest.approx(-1.0, abs=1e-10) + assert reals[1] == pytest.approx(1.0, abs=1e-10) + + def test_empty_system(self) -> None: + assert eigenvalues([]) == [] + + def test_complex_eigenvalues(self) -> None: + # [[0, 1], [-1, 0]] has eigenvalues ±j + eigs = eigenvalues([[0.0, 1.0], [-1.0, 0.0]]) + assert len(eigs) == 2 + for e in eigs: + assert abs(e.real) < 1e-10 + assert abs(abs(e.imag) - 1.0) < 1e-10 + + +class TestIsStable: + def test_stable_continuous(self) -> None: + assert is_stable([[-1.0, 0.0], [0.0, -2.0]]) + + def test_unstable_continuous(self) -> None: + assert not is_stable([[1.0, 0.0], [0.0, -1.0]]) + + def test_marginally_stable_is_not_asymptotically_stable(self) -> None: + # Pure oscillator: eigenvalues on imaginary axis + assert not is_stable([[0.0, 1.0], [-1.0, 0.0]]) + + def test_stable_discrete(self) -> None: + # |λ| < 1 for discrete stability + assert is_stable([[0.5, 0.0], [0.0, 0.3]], continuous=False) + + def test_unstable_discrete(self) -> None: + assert not is_stable([[1.1, 0.0], [0.0, 0.5]], continuous=False) + + def test_empty_system(self) -> None: + assert is_stable([]) + + +class TestIsMarginallyStable: + def test_oscillator(self) -> None: + assert is_marginally_stable([[0.0, 1.0], [-1.0, 0.0]]) + + def test_stable_is_not_marginal(self) -> None: + assert not is_marginally_stable([[-1.0, 0.0], [0.0, -2.0]]) + + +# --------------------------------------------------------------------------- +# Frequency response +# --------------------------------------------------------------------------- + + +class TestFrequencyResponse: + def test_first_order_bode(self) -> None: + """1/(s+1): -3dB at w=1, -45deg at w=1.""" + omega, mag_db, phase_deg = frequency_response( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + omega=[1.0], + ) + + assert len(omega) == 1 + assert mag_db[0] == pytest.approx(-3.01, abs=0.1) + assert phase_deg[0] == pytest.approx(-45.0, abs=2.0) + + def test_auto_frequency_range(self) -> None: + omega, mag_db, phase_deg = frequency_response( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + ) + assert len(omega) == 500 # default n_points + + +class TestGainMargin: + def test_known_system(self) -> None: + """1/((s+1)(s+2)(s+3)): known gain margin.""" + # H(s) = 1/(s^3 + 6s^2 + 11s + 6) + num = [1.0] + den = [1.0, 6.0, 11.0, 6.0] + gm_db, gm_freq = gain_margin(num, den) + assert gm_db > 0 # System is stable, should have positive margin + + +class TestPhaseMargin: + def test_known_system(self) -> None: + """1/(s+1): infinite phase margin (gain never reaches 0dB at + a frequency where phase is problematic).""" + num = [1.0] + den = [1.0, 1.0] + pm_deg, pm_freq = phase_margin(num, den) + # For 1/(s+1), gain is always < 0dB, so phase margin is infinite + assert pm_deg == float("inf") + + +# --------------------------------------------------------------------------- +# Discretization +# --------------------------------------------------------------------------- + + +class TestDiscretize: + def test_tustin_preserves_stability(self) -> None: + """Stable continuous system → stable discrete system.""" + Ad, Bd, Cd, Dd = discretize( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + dt=0.1, method="tustin", + ) + assert is_stable(Ad, continuous=False) + + def test_zoh(self) -> None: + """ZOH discretization of 1/(s+1) with dt=0.1.""" + Ad, Bd, Cd, Dd = discretize( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + dt=0.1, method="zoh", + ) + # A_d = e^{A*dt} = e^{-0.1} ≈ 0.9048 + assert Ad[0][0] == pytest.approx(0.9048, abs=0.01) + assert is_stable(Ad, continuous=False) + + def test_euler(self) -> None: + Ad, Bd, Cd, Dd = discretize( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + dt=0.1, method="euler", + ) + # Forward Euler: Ad = I + dt*A = [[0.9]] + assert Ad[0][0] == pytest.approx(0.9, abs=0.01) + + def test_invalid_method(self) -> None: + with pytest.raises(ValueError, match="Unknown discretization method"): + discretize( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + dt=0.1, method="invalid", + ) + + +# --------------------------------------------------------------------------- +# Controller synthesis +# --------------------------------------------------------------------------- + + +class TestLQR: + def test_double_integrator(self) -> None: + """LQR on double integrator should produce stable closed-loop.""" + A = [[0.0, 1.0], [0.0, 0.0]] + B = [[0.0], [1.0]] + Q = [[1.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + + K, P, E = lqr(A, B, Q, R) + + # K should be 1x2 + assert len(K) == 1 + assert len(K[0]) == 2 + + # Closed-loop should be stable + assert all(e.real < 0 for e in E) + + def test_scalar_system(self) -> None: + """Simple scalar system: dx/dt = -x + u, Q=1, R=1.""" + K, P, E = lqr( + A=[[-1.0]], B=[[1.0]], Q=[[1.0]], R=[[1.0]], + ) + assert len(K) == 1 + assert len(K[0]) == 1 + assert E[0].real < 0 + + +class TestDLQR: + def test_discrete_double_integrator(self) -> None: + """Discrete LQR on discretized double integrator.""" + # First discretize + Ad, Bd, Cd, Dd = discretize( + A=[[0.0, 1.0], [0.0, 0.0]], + B=[[0.0], [1.0]], + C=[[1.0, 0.0]], + D=[[0.0]], + dt=0.1, method="zoh", + ) + + Q = [[1.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + + K, P, E = dlqr(Ad, Bd, Q, R) + + # Closed-loop eigenvalues inside unit circle + assert all(abs(e) < 1.0 for e in E) + + +class TestKalman: + def test_observer_gain(self) -> None: + """Kalman filter for simple observable system.""" + A = [[-1.0, 0.0], [0.0, -2.0]] + C = [[1.0, 0.0]] + Q_proc = [[1.0, 0.0], [0.0, 1.0]] + R_meas = [[0.1]] + + L, P = kalman(A, C, Q_proc, R_meas) + + # L should be 2x1 (n x p) + assert len(L) == 2 + assert len(L[0]) == 1 + + # Observer poles = eigenvalues of A - LC should be stable + import numpy as np + + A_obs = np.array(A) - np.array(L) @ np.array(C) + eigs = np.linalg.eigvals(A_obs) + assert all(e.real < 0 for e in eigs) diff --git a/packages/gds-analysis/tests/test_response.py b/packages/gds-analysis/tests/test_response.py new file mode 100644 index 0000000..0c56104 --- /dev/null +++ b/packages/gds-analysis/tests/test_response.py @@ -0,0 +1,152 @@ +"""Tests for step/impulse response computation and metrics.""" + +import math + +import pytest + +from gds_analysis.response import StepMetrics, step_response_metrics + + +class TestStepResponseMetrics: + def test_first_order_no_overshoot(self) -> None: + """First-order system 1/(s+1): no overshoot, known time constants.""" + # Generate 1 - e^{-t} response + n = 1000 + t_max = 10.0 + times = [i * t_max / (n - 1) for i in range(n)] + values = [1.0 - math.exp(-t) for t in times] + + metrics = step_response_metrics(times, values, setpoint=1.0) + + # Rise time (10% to 90%): for 1-e^{-t}, t_10 = -ln(0.9) ≈ 0.105, + # t_90 = -ln(0.1) ≈ 2.303, rise time ≈ 2.197 + assert metrics.rise_time == pytest.approx(2.197, abs=0.05) + + # No overshoot + assert metrics.overshoot_pct == pytest.approx(0.0, abs=0.1) + + # Settling time (2% band): 1 - e^{-t} within 0.02 of 1.0 + # e^{-t} < 0.02 → t > -ln(0.02) ≈ 3.912 + assert metrics.settling_time == pytest.approx(3.91, abs=0.15) + + # Steady-state + assert metrics.steady_state_value == pytest.approx(1.0, abs=0.01) + assert metrics.steady_state_error == pytest.approx(0.0, abs=0.01) + + def test_underdamped_overshoot(self) -> None: + """Second-order underdamped: known overshoot formula.""" + zeta = 0.3 + wn = 2.0 + wd = wn * math.sqrt(1 - zeta**2) + + n = 2000 + t_max = 10.0 + times = [i * t_max / (n - 1) for i in range(n)] + values = [ + 1.0 - math.exp(-zeta * wn * t) / math.sqrt(1 - zeta**2) + * math.sin(wd * t + math.acos(zeta)) + for t in times + ] + + metrics = step_response_metrics(times, values, setpoint=1.0) + + # Analytical overshoot: exp(-pi*zeta/sqrt(1-zeta^2)) * 100 + expected_os = math.exp(-math.pi * zeta / math.sqrt(1 - zeta**2)) * 100 + assert metrics.overshoot_pct == pytest.approx(expected_os, abs=2.0) + + def test_steady_state_error(self) -> None: + """System that settles to 0.95 instead of 1.0.""" + n = 500 + times = [i * 10.0 / (n - 1) for i in range(n)] + values = [0.95 * (1.0 - math.exp(-t)) for t in times] + + metrics = step_response_metrics(times, values, setpoint=1.0) + assert metrics.steady_state_error == pytest.approx(0.05, abs=0.02) + + def test_too_few_points(self) -> None: + with pytest.raises(ValueError, match="at least 2"): + step_response_metrics([0.0], [0.0]) + + def test_mismatched_lengths(self) -> None: + with pytest.raises(ValueError, match="same length"): + step_response_metrics([0.0, 1.0], [0.0]) + + def test_constant_response(self) -> None: + """Constant response — degenerate but shouldn't crash.""" + times = [0.0, 1.0, 2.0, 3.0, 4.0] + values = [1.0, 1.0, 1.0, 1.0, 1.0] + metrics = step_response_metrics(times, values, setpoint=1.0) + assert metrics.rise_time == 0.0 + assert metrics.overshoot_pct == 0.0 + + def test_dataclass_fields(self) -> None: + """Verify StepMetrics has all required fields.""" + m = StepMetrics( + rise_time=1.0, + settling_time=2.0, + overshoot_pct=10.0, + peak_time=1.5, + peak_value=1.1, + steady_state_value=1.0, + steady_state_error=0.0, + ) + assert m.rise_time == 1.0 + assert m.settling_time == 2.0 + assert m.overshoot_pct == 10.0 + + +class TestStepResponse: + """Tests for step_response() — requires scipy.""" + + def test_first_order_final_value(self) -> None: + """1/(s+1) step response should approach 1.0.""" + from gds_analysis.response import step_response + + times, outputs = step_response( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + t_span=(0.0, 10.0), n_points=500, + ) + + assert len(times) == 500 + assert len(outputs) == 1 + # Final value should be close to 1.0 + assert outputs[0][-1] == pytest.approx(1.0, abs=0.01) + # Initial value should be close to 0.0 + assert outputs[0][0] == pytest.approx(0.0, abs=0.01) + + def test_double_integrator_ramp(self) -> None: + """1/s^2 step response is a parabola (t^2/2).""" + from gds_analysis.response import step_response + + times, outputs = step_response( + A=[[0.0, 1.0], [0.0, 0.0]], + B=[[0.0], [1.0]], + C=[[1.0, 0.0]], + D=[[0.0]], + t_span=(0.0, 5.0), n_points=500, + ) + + # At t=1, output ≈ 0.5; at t=2, output ≈ 2.0 + idx_t1 = 100 # t ≈ 1.0 + idx_t2 = 200 # t ≈ 2.0 + t1 = times[idx_t1] + t2 = times[idx_t2] + assert outputs[0][idx_t1] == pytest.approx(t1**2 / 2, abs=0.05) + assert outputs[0][idx_t2] == pytest.approx(t2**2 / 2, abs=0.1) + + +class TestImpulseResponse: + def test_first_order_decay(self) -> None: + """1/(s+1) impulse response = e^{-t}.""" + from gds_analysis.response import impulse_response + + times, outputs = impulse_response( + A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + t_span=(0.0, 10.0), n_points=500, + ) + + assert len(times) == 500 + # Initial value should be close to 1.0 + assert outputs[0][0] == pytest.approx(1.0, abs=0.1) + # Should decay toward 0 + assert outputs[0][-1] == pytest.approx(0.0, abs=0.01) diff --git a/packages/gds-domains/gds_domains/symbolic/__init__.py b/packages/gds-domains/gds_domains/symbolic/__init__.py index ea2a131..b9eb879 100644 --- a/packages/gds-domains/gds_domains/symbolic/__init__.py +++ b/packages/gds-domains/gds_domains/symbolic/__init__.py @@ -11,6 +11,20 @@ ) from gds_domains.symbolic.linearize import LinearizedSystem from gds_domains.symbolic.model import SymbolicControlModel +from gds_domains.symbolic.transfer import ( + TransferFunction, + TransferFunctionMatrix, + characteristic_polynomial, + controllability_matrix, + is_controllable, + is_minimum_phase, + is_observable, + observability_matrix, + poles, + sensitivity, + ss_to_tf, + zeros, +) __all__ = [ "HamiltonianSpec", @@ -20,7 +34,19 @@ "StateEquation", "SymbolicControlModel", "SymbolicError", + "TransferFunction", + "TransferFunctionMatrix", + "characteristic_polynomial", + "controllability_matrix", "derive_from_model", "derive_hamiltonian", + "is_controllable", + "is_minimum_phase", + "is_observable", + "observability_matrix", + "poles", + "sensitivity", + "ss_to_tf", "verify_conservation", + "zeros", ] diff --git a/packages/gds-domains/gds_domains/symbolic/delay.py b/packages/gds-domains/gds_domains/symbolic/delay.py new file mode 100644 index 0000000..7d77acf --- /dev/null +++ b/packages/gds-domains/gds_domains/symbolic/delay.py @@ -0,0 +1,120 @@ +"""Padé approximation for time delay modeling. + +Provides rational polynomial approximations of pure time delays e^{-sτ}, +enabling frequency-domain analysis of systems with transport delays. + +Uses SymPy only — no scipy or numpy dependency. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from gds_domains.symbolic._compat import require_sympy + +if TYPE_CHECKING: + from gds_domains.symbolic.transfer import TransferFunction + + +def pade_approximation(delay: float, order: int = 3) -> TransferFunction: + """Compute an (N, N) Padé approximation of e^{-sτ}. + + The Padé approximant is an all-pass rational function that matches + the Taylor expansion of e^{-sτ} up to order 2N+1. It has N zeros + in the right half-plane and N poles in the left half-plane. + + Parameters + ---------- + delay : float + Time delay τ in seconds. Must be > 0. + order : int + Approximation order N (1..8 typical). Higher order = more accurate + but adds RHP zeros that complicate root locus analysis. + + Returns + ------- + TransferFunction + Rational approximation with ``num`` and ``den`` coefficient lists. + + Raises + ------ + ValueError + If delay <= 0 or order < 1. + """ + from gds_domains.symbolic.transfer import TransferFunction + + if delay <= 0: + raise ValueError(f"delay must be > 0, got {delay}") + if order < 1: + raise ValueError(f"order must be >= 1, got {order}") + + require_sympy() + import sympy + + s = sympy.Symbol("s") + tau = sympy.Rational(delay).limit_denominator(10**12) + + # Padé coefficients using the closed form: + # P_n(s) = sum_{k=0}^{n} (-1)^k * C(n,k) * (s*tau)^k * (2n-k)! * n! + # / ((2n)! * k!) + # Q_n(s) = sum_{k=0}^{n} C(n,k) * (s*tau)^k * (2n-k)! * n! + # / ((2n)! * k!) + # + # Simplified: coefficient of (s*tau)^k in numerator and denominator: + # a_k = (2n - k)! * n! / ((2n)! * k! * (n - k)!) + # Numerator has alternating signs. + + n = order + num_poly = sympy.Integer(0) + den_poly = sympy.Integer(0) + + for k in range(n + 1): + coeff = ( + sympy.factorial(2 * n - k) + * sympy.factorial(n) + / (sympy.factorial(2 * n) * sympy.factorial(k) * sympy.factorial(n - k)) + ) + term = coeff * (s * tau) ** k + den_poly += term + num_poly += (-1) ** k * term + + num_sympy = sympy.Poly(num_poly, s) + den_sympy = sympy.Poly(den_poly, s) + + return TransferFunction( + num=[float(c) for c in num_sympy.all_coeffs()], + den=[float(c) for c in den_sympy.all_coeffs()], + ) + + +def delay_system( + tf: TransferFunction, + delay: float, + order: int = 3, +) -> TransferFunction: + """Cascade a transfer function with a Padé-approximated delay. + + Computes H_delayed(s) = H(s) * Padé(s, τ) by multiplying + the numerator and denominator polynomials. + + Parameters + ---------- + tf : TransferFunction + Original transfer function. + delay : float + Time delay τ in seconds. + order : int + Padé approximation order. + + Returns + ------- + TransferFunction + H(s) * Padé(s, τ) as a single transfer function. + """ + from gds_domains.symbolic.transfer import _tf_multiply + + if delay <= 0: + raise ValueError(f"delay must be > 0, got {delay}") + + pade_tf = pade_approximation(delay, order) + return _tf_multiply(tf, pade_tf) diff --git a/packages/gds-domains/gds_domains/symbolic/transfer.py b/packages/gds-domains/gds_domains/symbolic/transfer.py new file mode 100644 index 0000000..62849c6 --- /dev/null +++ b/packages/gds-domains/gds_domains/symbolic/transfer.py @@ -0,0 +1,437 @@ +"""Transfer function representation and classical control analysis. + +Provides symbolic transfer function computation from state-space (A, B, C, D) +matrices, pole/zero analysis, controllability/observability tests, and +sensitivity functions (Gang of Six). + +Uses SymPy only — no scipy or numpy dependency. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from gds_domains.symbolic._compat import require_sympy + +if TYPE_CHECKING: + from gds_domains.symbolic.linearize import LinearizedSystem + + +@dataclass(frozen=True) +class TransferFunction: + """SISO transfer function as coefficient lists (descending powers of s). + + ``num = [b_n, b_{n-1}, ..., b_0]`` and ``den = [a_m, a_{m-1}, ..., a_0]`` + represent H(s) = (b_n s^n + ... + b_0) / (a_m s^m + ... + a_0). + """ + + num: list[float] + den: list[float] + input_name: str = "" + output_name: str = "" + + +@dataclass(frozen=True) +class TransferFunctionMatrix: + """MIMO transfer function matrix (p outputs x m inputs). + + ``elements[i][j]`` is the SISO transfer function from input j to output i. + """ + + elements: list[list[TransferFunction]] + input_names: list[str] = field(default_factory=list) + output_names: list[str] = field(default_factory=list) + + +def ss_to_tf(ls: LinearizedSystem) -> TransferFunctionMatrix: + """Convert state-space (A, B, C, D) to a transfer function matrix. + + Computes H(s) = C(sI - A)^{-1}B + D using the adjugate/determinant + method for numerical stability with symbolic computation. + + Parameters + ---------- + ls : LinearizedSystem + State-space matrices with A, B, C, D as ``list[list[float]]``. + + Returns + ------- + TransferFunctionMatrix + p x m matrix of SISO transfer functions. + """ + require_sympy() + import sympy + + s = sympy.Symbol("s") + n = len(ls.A) + p = len(ls.C) # number of outputs + m = len(ls.B[0]) if n > 0 else len(ls.D[0]) if p > 0 else 0 + + A = sympy.Matrix(ls.A) + B = sympy.Matrix(ls.B) + C = sympy.Matrix(ls.C) + D = sympy.Matrix(ls.D) + + # sI - A + sI_A = s * sympy.eye(n) - A + + # Characteristic polynomial: det(sI - A) + char_poly = sI_A.det().as_poly(s) if n > 0 else sympy.Poly(1, s) + + # Adjugate of (sI - A) + adj = sI_A.adjugate() if n > 0 else sympy.Matrix() + + elements: list[list[TransferFunction]] = [] + for i in range(p): + row: list[TransferFunction] = [] + for j in range(m): + if n > 0: + # H_ij(s) = [C[i,:] @ adj @ B[:,j]] / det(sI-A) + D[i,j] + c_row = C.row(i) + b_col = B.col(j) + num_expr = (c_row * adj * b_col)[0, 0] + d_val = D[i, j] + + # H_ij = num_expr/char_poly + d_val + # = (num_expr + d_val * char_poly) / char_poly + num_poly = sympy.Poly(num_expr, s) + char_poly * d_val + den_poly = char_poly + else: + # No states — pure feedthrough + num_poly = sympy.Poly(D[i, j], s) + den_poly = sympy.Poly(1, s) + + # Extract coefficients (descending powers) + num_coeffs = [float(c) for c in num_poly.all_coeffs()] + den_coeffs = [float(c) for c in den_poly.all_coeffs()] + + row.append( + TransferFunction( + num=num_coeffs, + den=den_coeffs, + input_name=ls.input_names[j] if j < len(ls.input_names) else "", + output_name=( + ls.output_names[i] if i < len(ls.output_names) else "" + ), + ) + ) + elements.append(row) + + return TransferFunctionMatrix( + elements=elements, + input_names=list(ls.input_names), + output_names=list(ls.output_names), + ) + + +def characteristic_polynomial(ls: LinearizedSystem) -> list[float]: + """Compute the characteristic polynomial det(sI - A) as coefficients. + + Returns coefficients in descending powers of s: + ``[1, a_{n-1}, ..., a_0]`` for an n x n state matrix. + """ + require_sympy() + import sympy + + n = len(ls.A) + if n == 0: + return [1.0] + + s = sympy.Symbol("s") + A = sympy.Matrix(ls.A) + poly = (s * sympy.eye(n) - A).det().as_poly(s) + return [float(c) for c in poly.all_coeffs()] + + +def poles(tf: TransferFunction) -> list[complex]: + """Compute the poles (roots of denominator polynomial). + + Returns a list of complex numbers, one per pole (with multiplicity). + """ + require_sympy() + import sympy + + s = sympy.Symbol("s") + poly = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(tf.den))), s) + roots = sympy.roots(poly, multiple=True) + return [complex(r) for r in roots] + + +def zeros(tf: TransferFunction) -> list[complex]: + """Compute the zeros (roots of numerator polynomial). + + Returns a list of complex numbers, one per zero (with multiplicity). + """ + require_sympy() + import sympy + + s = sympy.Symbol("s") + poly = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(tf.num))), s) + roots = sympy.roots(poly, multiple=True) + return [complex(r) for r in roots] + + +def is_minimum_phase(tf: TransferFunction) -> bool: + """Check if all zeros have strictly negative real part. + + A minimum-phase system has no right-half-plane (RHP) zeros. + Returns True for systems with no zeros. + """ + z = zeros(tf) + if not z: + return True + return all(zi.real < 0 for zi in z) + + +def controllability_matrix(ls: LinearizedSystem) -> list[list[float]]: + """Compute the controllability matrix C_c = [B, AB, A^2B, ..., A^{n-1}B]. + + Returns an n x (n*m) matrix as ``list[list[float]]``. + """ + require_sympy() + import sympy + + n = len(ls.A) + if n == 0: + return [] + + A = sympy.Matrix(ls.A) + B = sympy.Matrix(ls.B) + + blocks = [B] + A_power = sympy.eye(n) + for _ in range(1, n): + A_power = A_power * A + blocks.append(A_power * B) + + # Concatenate horizontally + result = blocks[0] + for blk in blocks[1:]: + result = result.row_join(blk) + + return [[float(result[i, j]) for j in range(result.cols)] for i in range(n)] + + +def observability_matrix(ls: LinearizedSystem) -> list[list[float]]: + """Compute the observability matrix C_o = [C; CA; CA^2; ...; CA^{n-1}]. + + Returns an (n*p) x n matrix as ``list[list[float]]``. + """ + require_sympy() + import sympy + + n = len(ls.A) + if n == 0: + return [] + + A = sympy.Matrix(ls.A) + C = sympy.Matrix(ls.C) + + blocks = [C] + A_power = sympy.eye(n) + for _ in range(1, n): + A_power = A_power * A + blocks.append(C * A_power) + + # Concatenate vertically + result = blocks[0] + for blk in blocks[1:]: + result = result.col_join(blk) + + return [ + [float(result[i, j]) for j in range(result.cols)] for i in range(result.rows) + ] + + +def is_controllable(ls: LinearizedSystem) -> bool: + """Check if the system (A, B) is controllable. + + True if rank(C_c) == n (number of states). + """ + require_sympy() + import sympy + + n = len(ls.A) + if n == 0: + return True + + cm = controllability_matrix(ls) + mat = sympy.Matrix(cm) + return mat.rank() == n + + +def is_observable(ls: LinearizedSystem) -> bool: + """Check if the system (A, C) is observable. + + True if rank(C_o) == n (number of states). + """ + require_sympy() + import sympy + + n = len(ls.A) + if n == 0: + return True + + om = observability_matrix(ls) + mat = sympy.Matrix(om) + return mat.rank() == n + + +def _tf_multiply(a: TransferFunction, b: TransferFunction) -> TransferFunction: + """Multiply two transfer functions (convolve num and den).""" + require_sympy() + import sympy + + s = sympy.Symbol("s") + + a_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(a.num))), s) + a_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(a.den))), s) + b_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(b.num))), s) + b_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(b.den))), s) + + result_num = a_num * b_num + result_den = a_den * b_den + + return TransferFunction( + num=[float(c) for c in result_num.all_coeffs()], + den=[float(c) for c in result_den.all_coeffs()], + ) + + +def _tf_add(a: TransferFunction, b: TransferFunction) -> TransferFunction: + """Add two transfer functions.""" + require_sympy() + import sympy + + s = sympy.Symbol("s") + + a_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(a.num))), s) + a_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(a.den))), s) + b_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(b.num))), s) + b_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(b.den))), s) + + # a/a_d + b/b_d = (a*b_d + b*a_d) / (a_d * b_d) + result_num = a_num * b_den + b_num * a_den + result_den = a_den * b_den + + return TransferFunction( + num=[float(c) for c in result_num.all_coeffs()], + den=[float(c) for c in result_den.all_coeffs()], + ) + + +def _tf_feedback( + forward: TransferFunction, + feedback: TransferFunction | None = None, +) -> TransferFunction: + """Compute closed-loop TF: forward / (1 + forward * feedback). + + If feedback is None, uses unity feedback (feedback = 1). + """ + require_sympy() + import sympy + + s = sympy.Symbol("s") + + g_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(forward.num))), s) + g_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(forward.den))), s) + + if feedback is None: + # Unity feedback: G / (1 + G) = g_num / (g_den + g_num) + result_num = g_num * g_den # g_num * 1 (for general formula) + result_den = g_den * g_den + g_num * g_den # den*(den+num) + # Simplify: G/(1+G) = g_num / (g_den + g_num) + result_num = g_num + result_den = g_den + g_num + else: + h_num = sympy.Poly( + sum(c * s**i for i, c in enumerate(reversed(feedback.num))), s + ) + h_den = sympy.Poly( + sum(c * s**i for i, c in enumerate(reversed(feedback.den))), s + ) + # G/(1+GH) = (g_num * h_den) / (g_den * h_den + g_num * h_num) + result_num = g_num * h_den + result_den = g_den * h_den + g_num * h_num + + return TransferFunction( + num=[float(c) for c in result_num.all_coeffs()], + den=[float(c) for c in result_den.all_coeffs()], + ) + + +def sensitivity( + plant: TransferFunction, + controller: TransferFunction, +) -> dict[str, TransferFunction]: + """Compute the Gang of Six sensitivity functions. + + Given plant P and controller K with loop transfer function L = P * K: + + - **S** = 1 / (1 + L) — sensitivity + - **T** = L / (1 + L) — complementary sensitivity + - **CS** = K / (1 + L) — control sensitivity (K*S) + - **PS** = P / (1 + L) — load sensitivity (P*S) + - **KS** = K*S — noise → control (same as CS) + - **KPS** = K*P / (1 + L) — input disturbance sensitivity (T) + + Returns + ------- + dict with keys "S", "T", "CS", "PS", "KS", "KPS". + """ + require_sympy() + import sympy + + s = sympy.Symbol("s") + + # Build symbolic polynomials + p_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(plant.num))), s) + p_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(plant.den))), s) + k_num = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(controller.num))), s) + k_den = sympy.Poly(sum(c * s**i for i, c in enumerate(reversed(controller.den))), s) + + # L = P * K: L_num/L_den = (p_num * k_num) / (p_den * k_den) + l_num = p_num * k_num + l_den = p_den * k_den + + # 1 + L = (l_den + l_num) / l_den + one_plus_l_num = l_den + l_num # numerator of (1 + L) + + def _make_tf(num_poly: sympy.Poly, den_poly: sympy.Poly) -> TransferFunction: + return TransferFunction( + num=[float(c) for c in num_poly.all_coeffs()], + den=[float(c) for c in den_poly.all_coeffs()], + ) + + # S = 1/(1+L) = l_den / (l_den + l_num) + s_tf = _make_tf(l_den, one_plus_l_num) + + # T = L/(1+L) = l_num / (l_den + l_num) + t_tf = _make_tf(l_num, one_plus_l_num) + + # CS = K * S = K/(1+L) = (k_num * p_den) / (k_den * (l_den + l_num)) + # Actually: K/(1+L) = (k_num/k_den) * (l_den / (l_den + l_num)) + # = (k_num * l_den) / (k_den * (l_den + l_num)) + # Since l_den = p_den * k_den: + # = (k_num * p_den * k_den) / (k_den * (l_den + l_num)) + # = (k_num * p_den) / (l_den + l_num) + cs_tf = _make_tf(k_num * p_den, one_plus_l_num) + + # PS = P * S = P/(1+L) = (p_num * k_den) / (l_den + l_num) + ps_tf = _make_tf(p_num * k_den, one_plus_l_num) + + # KS = K * S (same as CS for standard feedback) + ks_tf = cs_tf + + # KPS = K * P * S = L/(1+L) = T + kps_tf = t_tf + + return { + "S": s_tf, + "T": t_tf, + "CS": cs_tf, + "PS": ps_tf, + "KS": ks_tf, + "KPS": kps_tf, + } diff --git a/packages/gds-domains/tests/test_delay.py b/packages/gds-domains/tests/test_delay.py new file mode 100644 index 0000000..3ab0b5b --- /dev/null +++ b/packages/gds-domains/tests/test_delay.py @@ -0,0 +1,99 @@ +"""Tests for Padé approximation of time delays.""" + +import math + +import pytest + +from gds_domains.symbolic.delay import delay_system, pade_approximation +from gds_domains.symbolic.transfer import TransferFunction, poles, zeros + + +class TestPadeApproximation: + def test_order_1_coefficients(self) -> None: + """Order-1 Padé of 1s delay: (−s/2 + 1)/(s/2 + 1).""" + tf = pade_approximation(1.0, order=1) + # num and den should be order 1 polynomials + assert len(tf.num) == 2 + assert len(tf.den) == 2 + + def test_allpass_property(self) -> None: + """Padé approximation should be approximately all-pass: |H(jw)| ≈ 1.""" + tf = pade_approximation(0.5, order=3) + + # Evaluate |H(jw)| at several frequencies using polynomial evaluation + for omega in [0.1, 1.0, 5.0, 10.0]: + jw = complex(0, omega) + num_val = sum( + c * jw**i for i, c in enumerate(reversed(tf.num)) + ) + den_val = sum( + c * jw**i for i, c in enumerate(reversed(tf.den)) + ) + mag = abs(num_val / den_val) + assert mag == pytest.approx(1.0, abs=0.01), ( + f"|H(j{omega})| = {mag}, expected ≈ 1.0" + ) + + def test_rhp_zeros(self) -> None: + """Order-N Padé should have N zeros in the right half-plane.""" + for order in [1, 2, 3]: + tf = pade_approximation(0.5, order=order) + z = zeros(tf) + rhp_zeros = [zi for zi in z if zi.real > 0] + assert len(rhp_zeros) == order + + def test_lhp_poles(self) -> None: + """Order-N Padé should have N poles in the left half-plane.""" + for order in [1, 2, 3]: + tf = pade_approximation(0.5, order=order) + p = poles(tf) + lhp_poles = [pi for pi in p if pi.real < 0] + assert len(lhp_poles) == order + + def test_phase_accuracy_low_frequency(self) -> None: + """At low frequencies, Padé phase should match -w*tau.""" + tau = 0.5 + tf = pade_approximation(tau, order=3) + + for omega in [0.1, 0.5, 1.0]: + jw = complex(0, omega) + num_val = sum( + c * jw**i for i, c in enumerate(reversed(tf.num)) + ) + den_val = sum( + c * jw**i for i, c in enumerate(reversed(tf.den)) + ) + H = num_val / den_val + pade_phase = math.atan2(H.imag, H.real) + exact_phase = -omega * tau + + assert pade_phase == pytest.approx(exact_phase, abs=0.05), ( + f"At w={omega}: Padé phase={pade_phase:.4f}, " + f"exact={exact_phase:.4f}" + ) + + def test_invalid_delay(self) -> None: + with pytest.raises(ValueError, match="delay must be > 0"): + pade_approximation(0.0) + with pytest.raises(ValueError, match="delay must be > 0"): + pade_approximation(-1.0) + + def test_invalid_order(self) -> None: + with pytest.raises(ValueError, match="order must be >= 1"): + pade_approximation(1.0, order=0) + + +class TestDelaySystem: + def test_pole_count(self) -> None: + """Delayed system should have original poles + Padé poles.""" + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) # 1/(s+1) + delayed = delay_system(plant, delay=0.5, order=2) + + p = poles(delayed) + # Original: 1 pole. Padé order 2: 2 poles. Total: 3. + assert len(p) == 3 + + def test_invalid_delay(self) -> None: + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) + with pytest.raises(ValueError, match="delay must be > 0"): + delay_system(plant, delay=0.0) diff --git a/packages/gds-domains/tests/test_transfer.py b/packages/gds-domains/tests/test_transfer.py new file mode 100644 index 0000000..63dd8d2 --- /dev/null +++ b/packages/gds-domains/tests/test_transfer.py @@ -0,0 +1,302 @@ +"""Tests for symbolic transfer function analysis.""" + +import pytest + +from gds_domains.symbolic.linearize import LinearizedSystem +from gds_domains.symbolic.transfer import ( + TransferFunction, + characteristic_polynomial, + controllability_matrix, + is_controllable, + is_minimum_phase, + is_observable, + observability_matrix, + poles, + sensitivity, + ss_to_tf, + zeros, +) + + +# --------------------------------------------------------------------------- +# Fixtures: canonical test systems +# --------------------------------------------------------------------------- + + +@pytest.fixture() +def first_order_system() -> LinearizedSystem: + """1/(s+1): A=[-1], B=[1], C=[1], D=[0].""" + return LinearizedSystem( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + x0=[0.0], + u0=[0.0], + state_names=["x"], + input_names=["u"], + output_names=["y"], + ) + + +@pytest.fixture() +def double_integrator() -> LinearizedSystem: + """1/s^2: two poles at origin.""" + return LinearizedSystem( + A=[[0.0, 1.0], [0.0, 0.0]], + B=[[0.0], [1.0]], + C=[[1.0, 0.0]], + D=[[0.0]], + x0=[0.0, 0.0], + u0=[0.0], + state_names=["position", "velocity"], + input_names=["force"], + output_names=["position_out"], + ) + + +@pytest.fixture() +def second_order_underdamped() -> LinearizedSystem: + """wn=1, zeta=0.5: poles at -0.5 +/- j*sqrt(3)/2.""" + return LinearizedSystem( + A=[[0.0, 1.0], [-1.0, -1.0]], + B=[[0.0], [1.0]], + C=[[1.0, 0.0]], + D=[[0.0]], + x0=[0.0, 0.0], + u0=[0.0], + state_names=["x", "v"], + input_names=["u"], + output_names=["y"], + ) + + +@pytest.fixture() +def non_minimum_phase() -> LinearizedSystem: + """System with a RHP zero: H(s) = (s-1)/((s+1)(s+2)).""" + return LinearizedSystem( + A=[[0.0, 1.0], [-2.0, -3.0]], + B=[[0.0], [1.0]], + C=[[-1.0, 1.0]], + D=[[0.0]], + x0=[0.0, 0.0], + u0=[0.0], + state_names=["x1", "x2"], + input_names=["u"], + output_names=["y"], + ) + + +@pytest.fixture() +def uncontrollable_system() -> LinearizedSystem: + """System where B doesn't span the state space.""" + return LinearizedSystem( + A=[[-1.0, 0.0], [0.0, -2.0]], + B=[[1.0], [0.0]], + C=[[1.0, 1.0]], + D=[[0.0]], + x0=[0.0, 0.0], + u0=[0.0], + state_names=["x1", "x2"], + input_names=["u"], + output_names=["y"], + ) + + +# --------------------------------------------------------------------------- +# Transfer function conversion +# --------------------------------------------------------------------------- + + +class TestSsToTf: + def test_first_order(self, first_order_system: LinearizedSystem) -> None: + tf_mat = ss_to_tf(first_order_system) + assert len(tf_mat.elements) == 1 + assert len(tf_mat.elements[0]) == 1 + tf = tf_mat.elements[0][0] + + # H(s) = 1/(s+1) → num=[1], den=[1, 1] + assert len(tf.num) == 1 + assert len(tf.den) == 2 + assert tf.num[0] == pytest.approx(1.0, abs=1e-10) + assert tf.den[0] == pytest.approx(1.0, abs=1e-10) + assert tf.den[1] == pytest.approx(1.0, abs=1e-10) + + def test_double_integrator(self, double_integrator: LinearizedSystem) -> None: + tf_mat = ss_to_tf(double_integrator) + tf = tf_mat.elements[0][0] + + # H(s) = 1/s^2 → num=[1], den=[1, 0, 0] + assert tf.num[-1] == pytest.approx(1.0, abs=1e-10) + assert tf.den[0] == pytest.approx(1.0, abs=1e-10) + assert tf.den[1] == pytest.approx(0.0, abs=1e-10) + assert tf.den[2] == pytest.approx(0.0, abs=1e-10) + + def test_preserves_names(self, first_order_system: LinearizedSystem) -> None: + tf_mat = ss_to_tf(first_order_system) + assert tf_mat.input_names == ["u"] + assert tf_mat.output_names == ["y"] + assert tf_mat.elements[0][0].input_name == "u" + assert tf_mat.elements[0][0].output_name == "y" + + +# --------------------------------------------------------------------------- +# Characteristic polynomial +# --------------------------------------------------------------------------- + + +class TestCharacteristicPolynomial: + def test_first_order(self, first_order_system: LinearizedSystem) -> None: + cp = characteristic_polynomial(first_order_system) + # det(sI - [-1]) = s + 1 + assert len(cp) == 2 + assert cp[0] == pytest.approx(1.0, abs=1e-10) + assert cp[1] == pytest.approx(1.0, abs=1e-10) + + def test_double_integrator(self, double_integrator: LinearizedSystem) -> None: + cp = characteristic_polynomial(double_integrator) + # det(sI - [[0,1],[0,0]]) = s^2 + assert len(cp) == 3 + assert cp[0] == pytest.approx(1.0, abs=1e-10) + assert cp[1] == pytest.approx(0.0, abs=1e-10) + assert cp[2] == pytest.approx(0.0, abs=1e-10) + + def test_empty_system(self) -> None: + ls = LinearizedSystem(A=[], B=[], C=[], D=[], x0=[], u0=[]) + cp = characteristic_polynomial(ls) + assert cp == [1.0] + + +# --------------------------------------------------------------------------- +# Poles and zeros +# --------------------------------------------------------------------------- + + +class TestPolesZeros: + def test_first_order_pole(self, first_order_system: LinearizedSystem) -> None: + tf_mat = ss_to_tf(first_order_system) + tf = tf_mat.elements[0][0] + p = poles(tf) + assert len(p) == 1 + assert p[0].real == pytest.approx(-1.0, abs=1e-6) + assert abs(p[0].imag) < 1e-6 + + def test_double_integrator_poles( + self, double_integrator: LinearizedSystem + ) -> None: + tf_mat = ss_to_tf(double_integrator) + tf = tf_mat.elements[0][0] + p = poles(tf) + assert len(p) == 2 + for pi in p: + assert abs(pi) < 1e-6 # Both at origin + + def test_underdamped_poles( + self, second_order_underdamped: LinearizedSystem + ) -> None: + tf_mat = ss_to_tf(second_order_underdamped) + tf = tf_mat.elements[0][0] + p = poles(tf) + assert len(p) == 2 + # Both should have negative real part + for pi in p: + assert pi.real < 0 + + def test_first_order_no_zeros(self, first_order_system: LinearizedSystem) -> None: + tf_mat = ss_to_tf(first_order_system) + tf = tf_mat.elements[0][0] + z = zeros(tf) + assert len(z) == 0 + + def test_minimum_phase(self, first_order_system: LinearizedSystem) -> None: + tf_mat = ss_to_tf(first_order_system) + assert is_minimum_phase(tf_mat.elements[0][0]) + + def test_non_minimum_phase(self, non_minimum_phase: LinearizedSystem) -> None: + tf_mat = ss_to_tf(non_minimum_phase) + tf = tf_mat.elements[0][0] + z = zeros(tf) + # Should have at least one RHP zero + assert any(zi.real > 0 for zi in z) + assert not is_minimum_phase(tf) + + +# --------------------------------------------------------------------------- +# Controllability / Observability +# --------------------------------------------------------------------------- + + +class TestControllabilityObservability: + def test_double_integrator_controllable( + self, double_integrator: LinearizedSystem + ) -> None: + assert is_controllable(double_integrator) + + def test_double_integrator_observable( + self, double_integrator: LinearizedSystem + ) -> None: + assert is_observable(double_integrator) + + def test_uncontrollable(self, uncontrollable_system: LinearizedSystem) -> None: + assert not is_controllable(uncontrollable_system) + + def test_controllability_matrix_shape( + self, double_integrator: LinearizedSystem + ) -> None: + cm = controllability_matrix(double_integrator) + assert len(cm) == 2 # n rows + assert len(cm[0]) == 2 # n*m columns (2 states, 1 input) + + def test_observability_matrix_shape( + self, double_integrator: LinearizedSystem + ) -> None: + om = observability_matrix(double_integrator) + assert len(om) == 2 # n*p rows (2 states, 1 output) + assert len(om[0]) == 2 # n columns + + def test_empty_system(self) -> None: + ls = LinearizedSystem(A=[], B=[], C=[], D=[], x0=[], u0=[]) + assert is_controllable(ls) + assert is_observable(ls) + + +# --------------------------------------------------------------------------- +# Sensitivity (Gang of Six) +# --------------------------------------------------------------------------- + + +class TestSensitivity: + def test_s_plus_t_equals_one(self) -> None: + """S + T = 1 is a fundamental identity.""" + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) # 1/(s+1) + controller = TransferFunction(num=[2.0], den=[1.0]) # K=2 + + gang = sensitivity(plant, controller) + + # Verify S + T = 1 by checking that num_S * den_T + num_T * den_S + # equals den_S * den_T (all should be the same denominator) + s_tf = gang["S"] + t_tf = gang["T"] + + # S and T should have the same denominator + assert len(s_tf.den) == len(t_tf.den) + for a, b in zip(s_tf.den, t_tf.den): + assert a == pytest.approx(b, abs=1e-10) + + # S_num + T_num should equal the denominator + # Pad shorter one + max_len = max(len(s_tf.num), len(t_tf.num)) + s_padded = [0.0] * (max_len - len(s_tf.num)) + s_tf.num + t_padded = [0.0] * (max_len - len(t_tf.num)) + t_tf.num + d_padded = [0.0] * (max_len - len(s_tf.den)) + s_tf.den + + summed = [a + b for a, b in zip(s_padded, t_padded)] + for a, b in zip(summed, d_padded): + assert a == pytest.approx(b, abs=1e-10) + + def test_returns_all_six(self) -> None: + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) + controller = TransferFunction(num=[1.0], den=[1.0]) + + gang = sensitivity(plant, controller) + assert set(gang.keys()) == {"S", "T", "CS", "PS", "KS", "KPS"} diff --git a/packages/gds-proof/gds_proof/analysis/lyapunov.py b/packages/gds-proof/gds_proof/analysis/lyapunov.py new file mode 100644 index 0000000..9cd52a8 --- /dev/null +++ b/packages/gds-proof/gds_proof/analysis/lyapunov.py @@ -0,0 +1,527 @@ +"""Lyapunov stability proofs and passivity certificates. + +Provides convenience constructors that translate control-theoretic stability +specifications into the invariant format the existing gds-proof engine handles. + +For continuous-time systems: + V(x) > 0 and dV/dt = (∂V/∂x) · f(x) < 0 + +For discrete-time systems: + V(x) > 0 and V(f(x)) - V(x) < 0 + +All proofs delegate to the existing five-strategy SymPy simplification +engine in ``gds_proof.analysis.symbolic``. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +import sympy +from sympy.parsing.sympy_parser import parse_expr + + +@dataclass(frozen=True) +class LyapunovResult: + """Result of a Lyapunov stability analysis.""" + + candidate: sympy.Expr + """V(x) expression.""" + + positive_definite: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + """Whether V(x) > 0 for x ≠ 0 was established.""" + + decreasing: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + """Whether dV/dt < 0 (continuous) or ΔV < 0 (discrete) was established.""" + + stable: bool + """True only if both positive_definite and decreasing are PROVED.""" + + dV_expr: sympy.Expr | None = None + """dV/dt or ΔV expression for inspection.""" + + details: list[str] = field(default_factory=list) + """Human-readable proof trace.""" + + +@dataclass(frozen=True) +class PassivityResult: + """Result of a passivity/dissipativity analysis.""" + + storage_function: sympy.Expr + """V(x) — storage function.""" + + supply_rate: sympy.Expr + """s(u, y) — supply rate.""" + + dissipation_proved: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + """Whether dV/dt ≤ s(u, y) was established.""" + + passive: bool + """True only if dissipation was PROVED.""" + + details: list[str] = field(default_factory=list) + + +def lyapunov_candidate( + V_expr: str | sympy.Expr, + state_transition: dict[str, str | sympy.Expr], + state_symbols: list[str], + *, + continuous: bool = True, +) -> LyapunovResult: + """Verify a Lyapunov candidate function. + + For continuous-time (``continuous=True``): + Checks V(x) > 0 and dV/dt = (∂V/∂x) · f(x) < 0. + + For discrete-time (``continuous=False``): + Checks V(x) > 0 and V(f(x)) - V(x) < 0. + + Parameters + ---------- + V_expr : str | sympy.Expr + Lyapunov candidate function. If str, parsed via SymPy. + state_transition : dict[str, str | sympy.Expr] + Maps state variable name to its dynamics. + Continuous: dx_i/dt = f_i(x). Discrete: x_i' = f_i(x). + state_symbols : list[str] + State variable names. + continuous : bool + Whether the system is continuous-time (True) or discrete-time (False). + + Returns + ------- + LyapunovResult + """ + details: list[str] = [] + + # Build symbol table + syms = {name: sympy.Symbol(name) for name in state_symbols} + + # Parse V + V = parse_expr(V_expr, local_dict=syms) if isinstance(V_expr, str) else V_expr + + # Parse state transition + f_exprs: dict[str, sympy.Expr] = {} + for name, expr in state_transition.items(): + if isinstance(expr, str): + f_exprs[name] = parse_expr(expr, local_dict=syms) + else: + f_exprs[name] = expr + + details.append(f"Candidate: V = {V}") + + # Check positive definiteness + # V(0) = 0 and V(x) > 0 for x ≠ 0 + zero_subs = {syms[n]: 0 for n in state_symbols} + V_at_zero = V.subs(zero_subs) + V_at_zero_simplified = sympy.simplify(V_at_zero) + + pd_status: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + if V_at_zero_simplified != 0: + pd_status = "FAILED" + details.append(f"V(0) = {V_at_zero_simplified} ≠ 0 → not positive definite") + else: + details.append("V(0) = 0 ✓") + # Try to verify V > 0 for x ≠ 0 using SymPy + # For polynomial V, check if it can be shown positive + pd_check = _check_positive_definite(V, [syms[n] for n in state_symbols]) + pd_status = pd_check + details.append(f"Positive definiteness: {pd_status}") + + # Check decrease condition + if continuous: + # dV/dt = sum_i (∂V/∂x_i) * f_i(x) + dV_dt = sympy.Integer(0) + for name in state_symbols: + partial = sympy.diff(V, syms[name]) + dV_dt += partial * f_exprs.get(name, sympy.Integer(0)) + dV_dt = sympy.expand(dV_dt) + details.append(f"dV/dt = {dV_dt}") + decrease_expr = dV_dt + else: + # ΔV = V(f(x)) - V(x) + subs_map = {syms[name]: f_exprs.get(name, syms[name]) for name in state_symbols} + V_next = V.subs(subs_map) + delta_V = sympy.expand(V_next - V) + details.append(f"ΔV = {delta_V}") + decrease_expr = delta_V + + # Check if decrease_expr < 0 + dec_status = _check_negative_definite( + decrease_expr, [syms[n] for n in state_symbols] + ) + details.append(f"Decrease condition: {dec_status}") + + stable = pd_status == "PROVED" and dec_status == "PROVED" + + return LyapunovResult( + candidate=V, + positive_definite=pd_status, + decreasing=dec_status, + stable=stable, + dV_expr=decrease_expr, + details=details, + ) + + +def quadratic_lyapunov( + P: list[list[float]], + A: list[list[float]], + state_symbols: list[str] | None = None, +) -> LyapunovResult: + """Verify V(x) = x'Px as a Lyapunov function for dx/dt = Ax. + + Checks: + 1. P is symmetric positive definite (all eigenvalues > 0). + 2. A'P + PA is negative definite (Lyapunov inequality). + + Parameters + ---------- + P : list[list[float]] + Candidate Lyapunov matrix (n x n, should be symmetric). + A : list[list[float]] + State matrix (n x n). + state_symbols : list[str] | None + State variable names. If None, uses x1, x2, ... + + Returns + ------- + LyapunovResult + """ + n = len(A) + if state_symbols is None: + state_symbols = [f"x{i + 1}" for i in range(n)] + + details: list[str] = [] + + P_mat = sympy.Matrix(P) + A_mat = sympy.Matrix(A) + + # Build V = x'Px symbolically + syms = [sympy.Symbol(name) for name in state_symbols] + x_vec = sympy.Matrix(syms) + V = (x_vec.T * P_mat * x_vec)[0, 0] + V = sympy.expand(V) + + details.append(f"V(x) = {V}") + + # Check P is symmetric positive definite via eigenvalues + eigenvals = P_mat.eigenvals() + eig_values = [] + for ev, mult in eigenvals.items(): + eig_values.extend([complex(ev)] * mult) + + pd_status: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + if all(e.real > 0 and abs(e.imag) < 1e-10 for e in eig_values): + pd_status = "PROVED" + details.append(f"P eigenvalues: {[float(e.real) for e in eig_values]} > 0 ✓") + else: + pd_status = "FAILED" + details.append(f"P eigenvalues: {eig_values} — not all positive") + + # Check A'P + PA is negative definite + Q_lyap = A_mat.T * P_mat + P_mat * A_mat + Q_eigenvals = Q_lyap.eigenvals() + q_eig_values = [] + for ev, mult in Q_eigenvals.items(): + q_eig_values.extend([complex(ev)] * mult) + + dec_status: Literal["PROVED", "FAILED", "INCONCLUSIVE"] + if all(e.real < 0 and abs(e.imag) < 1e-10 for e in q_eig_values): + dec_status = "PROVED" + details.append( + f"A'P + PA eigenvalues: {[float(e.real) for e in q_eig_values]} < 0 ✓" + ) + elif any(e.real > 0 for e in q_eig_values): + dec_status = "FAILED" + details.append(f"A'P + PA eigenvalues: {q_eig_values} — not all negative") + else: + dec_status = "INCONCLUSIVE" + details.append(f"A'P + PA eigenvalues: {q_eig_values} — inconclusive") + + # Compute dV/dt for reporting + dV_dt = (x_vec.T * Q_lyap * x_vec)[0, 0] + dV_dt = sympy.expand(dV_dt) + details.append(f"dV/dt = x'(A'P + PA)x = {dV_dt}") + + stable = pd_status == "PROVED" and dec_status == "PROVED" + + return LyapunovResult( + candidate=V, + positive_definite=pd_status, + decreasing=dec_status, + stable=stable, + dV_expr=dV_dt, + details=details, + ) + + +def find_quadratic_lyapunov( + A: list[list[float]], + Q: list[list[float]] | None = None, +) -> tuple[list[list[float]], LyapunovResult] | None: + """Attempt to find P satisfying A'P + PA = -Q (Lyapunov equation). + + If Q is None, uses Q = I (identity). + + Parameters + ---------- + A : list[list[float]] + State matrix. + Q : list[list[float]] | None + Desired negative-definiteness target. Default is identity. + + Returns + ------- + (P, result) if a valid P was found, None if the system is not stable. + """ + n = len(A) + if n == 0: + return None + + A_mat = sympy.Matrix(A) + + Q_mat = sympy.eye(n) if Q is None else sympy.Matrix(Q) + + # Solve A'P + PA = -Q by parameterizing P as symmetric + # and solving the resulting linear system + p_vars = {} + P_sym = sympy.zeros(n) + for i in range(n): + for j in range(i, n): + var = sympy.Symbol(f"p_{i}_{j}") + p_vars[(i, j)] = var + P_sym[i, j] = var + P_sym[j, i] = var + + # A'P + PA + Q = 0 + equation_matrix = A_mat.T * P_sym + P_sym * A_mat + Q_mat + + # Collect all equations + equations = [] + for i in range(n): + for j in range(i, n): + equations.append(equation_matrix[i, j]) + + variables = list(p_vars.values()) + + try: + solution = sympy.solve(equations, variables, dict=True) + except (NotImplementedError, ValueError): + return None + + if not solution: + return None + + sol = solution[0] + + # Reconstruct P + P_result = [[0.0] * n for _ in range(n)] + for (i, j), var in p_vars.items(): + val = float(sol.get(var, 0)) + P_result[i][j] = val + P_result[j][i] = val + + # Verify the result + result = quadratic_lyapunov(P_result, A) + + if result.stable: + return P_result, result + return None + + +def passivity_certificate( + V_expr: str | sympy.Expr, + supply_rate: str | sympy.Expr, + state_transition: dict[str, str | sympy.Expr], + output_map: dict[str, str | sympy.Expr], + state_symbols: list[str], + input_symbols: list[str], +) -> PassivityResult: + """Verify passivity: dV/dt ≤ s(u, y). + + For a passive system with supply rate s(u, y) = u'y: + dV/dt ≤ u'y (energy dissipation inequality) + + More general supply rates (e.g., s = γ²|u|² - |y|² for L2-gain) + can be specified. + + Parameters + ---------- + V_expr : str | sympy.Expr + Storage function V(x). + supply_rate : str | sympy.Expr + Supply rate s(u, y) as expression over input and output symbols. + state_transition : dict[str, str | sympy.Expr] + dx_i/dt = f_i(x, u). + output_map : dict[str, str | sympy.Expr] + y_i = h_i(x, u). + state_symbols : list[str] + State variable names. + input_symbols : list[str] + Input variable names. + """ + details: list[str] = [] + + # Build symbol table + all_sym_names = state_symbols + input_symbols + syms = {name: sympy.Symbol(name) for name in all_sym_names} + + # Add output symbols for supply rate parsing + out_syms = {name: sympy.Symbol(name) for name in output_map} + all_syms = {**syms, **out_syms} + + # Parse expressions + V = parse_expr(V_expr, local_dict=syms) if isinstance(V_expr, str) else V_expr + + if isinstance(supply_rate, str): + s_rate = parse_expr(supply_rate, local_dict=all_syms) + else: + s_rate = supply_rate + + f_exprs: dict[str, sympy.Expr] = {} + for name, expr in state_transition.items(): + if isinstance(expr, str): + f_exprs[name] = parse_expr(expr, local_dict=syms) + else: + f_exprs[name] = expr + + h_exprs: dict[str, sympy.Expr] = {} + for name, expr in output_map.items(): + if isinstance(expr, str): + h_exprs[name] = parse_expr(expr, local_dict=syms) + else: + h_exprs[name] = expr + + details.append(f"Storage function: V = {V}") + details.append(f"Supply rate: s = {s_rate}") + + # Compute dV/dt + dV_dt = sympy.Integer(0) + for name in state_symbols: + partial = sympy.diff(V, syms[name]) + dV_dt += partial * f_exprs.get(name, sympy.Integer(0)) + dV_dt = sympy.expand(dV_dt) + details.append(f"dV/dt = {dV_dt}") + + # Substitute output expressions into supply rate + out_subs = {out_syms[name]: h_exprs[name] for name in output_map} + s_rate_expanded = s_rate.subs(out_subs) + s_rate_expanded = sympy.expand(s_rate_expanded) + details.append(f"Supply rate (expanded): s = {s_rate_expanded}") + + # Check dV/dt - s(u, y) ≤ 0 + dissipation = sympy.expand(dV_dt - s_rate_expanded) + details.append(f"Dissipation: dV/dt - s = {dissipation}") + + # Try to prove dissipation ≤ 0 + result = _check_nonpositive(dissipation, list(syms.values())) + + details.append(f"Dissipation check: {result}") + + return PassivityResult( + storage_function=V, + supply_rate=s_rate, + dissipation_proved=result, + passive=result == "PROVED", + details=details, + ) + + +# --------------------------------------------------------------------------- +# Internal helpers +# --------------------------------------------------------------------------- + + +def _check_positive_definite( + expr: sympy.Expr, + symbols: list[sympy.Symbol], +) -> Literal["PROVED", "FAILED", "INCONCLUSIVE"]: + """Attempt to verify that expr > 0 for all non-zero symbol values. + + Uses Hessian check for quadratic forms and SymPy simplification. + """ + # Check if it's a quadratic form x'Hx + # Compute the Hessian matrix + n = len(symbols) + if n == 0: + return "INCONCLUSIVE" + + hessian = sympy.zeros(n) + for i in range(n): + for j in range(n): + hessian[i, j] = sympy.diff(expr, symbols[i], symbols[j]) + + # Check if Hessian is constant (quadratic form) + is_constant_hessian = all( + hessian[i, j].free_symbols == set() for i in range(n) for j in range(n) + ) + + if is_constant_hessian: + # For V = (1/2) x'Hx, H must be positive definite + H = hessian + eigenvals = H.eigenvals() + eig_values = [] + for ev, mult in eigenvals.items(): + eig_values.extend([complex(ev)] * mult) + + if all(e.real > 0 and abs(e.imag) < 1e-10 for e in eig_values): + return "PROVED" + if any(e.real < 0 for e in eig_values): + return "FAILED" + + return "INCONCLUSIVE" + + +def _check_negative_definite( + expr: sympy.Expr, + symbols: list[sympy.Symbol], +) -> Literal["PROVED", "FAILED", "INCONCLUSIVE"]: + """Attempt to verify that expr < 0 for all non-zero symbol values.""" + result = _check_positive_definite(-expr, symbols) + return result + + +def _check_nonpositive( + expr: sympy.Expr, + symbols: list[sympy.Symbol], +) -> Literal["PROVED", "FAILED", "INCONCLUSIVE"]: + """Attempt to verify that expr <= 0 for all symbol values.""" + # First try negative definiteness + result = _check_negative_definite(expr, symbols) + if result == "PROVED": + return "PROVED" + + # Check if expr simplifies to 0 + simplified = sympy.simplify(expr) + if simplified == 0: + return "PROVED" + + # Check negative semi-definiteness via Hessian eigenvalues + n = len(symbols) + if n > 0: + hessian = sympy.zeros(n) + for i in range(n): + for j in range(n): + hessian[i, j] = sympy.diff(expr, symbols[i], symbols[j]) + + is_constant_hessian = all( + hessian[i, j].free_symbols == set() for i in range(n) for j in range(n) + ) + + if is_constant_hessian: + eigenvals = hessian.eigenvals() + eig_values = [] + for ev, mult in eigenvals.items(): + eig_values.extend([complex(ev)] * mult) + + # Negative semi-definite: all eigenvalues <= 0 + if all(e.real <= 1e-10 and abs(e.imag) < 1e-10 for e in eig_values): + return "PROVED" + if any(e.real > 1e-10 for e in eig_values): + return "FAILED" + + return result diff --git a/packages/gds-proof/tests/test_lyapunov.py b/packages/gds-proof/tests/test_lyapunov.py new file mode 100644 index 0000000..929e114 --- /dev/null +++ b/packages/gds-proof/tests/test_lyapunov.py @@ -0,0 +1,175 @@ +"""Tests for Lyapunov stability proofs and passivity certificates.""" + +import pytest + +from gds_proof.analysis.lyapunov import ( + LyapunovResult, + PassivityResult, + find_quadratic_lyapunov, + lyapunov_candidate, + passivity_certificate, + quadratic_lyapunov, +) + + +class TestLyapunovCandidate: + def test_stable_linear_continuous(self) -> None: + """dx/dt = -x, V = x^2: should prove stable.""" + result = lyapunov_candidate( + V_expr="x**2", + state_transition={"x": "-x"}, + state_symbols=["x"], + continuous=True, + ) + assert result.positive_definite == "PROVED" + assert result.decreasing == "PROVED" + assert result.stable is True + + def test_unstable_linear_continuous(self) -> None: + """dx/dt = x, V = x^2: dV/dt = 2x^2 > 0, should fail decrease.""" + result = lyapunov_candidate( + V_expr="x**2", + state_transition={"x": "x"}, + state_symbols=["x"], + continuous=True, + ) + assert result.positive_definite == "PROVED" + assert result.decreasing == "FAILED" + assert result.stable is False + + def test_2d_stable(self) -> None: + """dx/dt = [-x, -2y], V = x^2 + y^2.""" + result = lyapunov_candidate( + V_expr="x**2 + y**2", + state_transition={"x": "-x", "y": "-2*y"}, + state_symbols=["x", "y"], + continuous=True, + ) + assert result.positive_definite == "PROVED" + assert result.decreasing == "PROVED" + assert result.stable is True + + def test_discrete_stable(self) -> None: + """x' = 0.5*x, V = x^2: V(f(x)) - V(x) = 0.25x^2 - x^2 = -0.75x^2 < 0.""" + result = lyapunov_candidate( + V_expr="x**2", + state_transition={"x": "0.5*x"}, + state_symbols=["x"], + continuous=False, + ) + assert result.positive_definite == "PROVED" + assert result.decreasing == "PROVED" + assert result.stable is True + + def test_discrete_unstable(self) -> None: + """x' = 2*x, V = x^2: V(f(x)) = 4x^2 > x^2.""" + result = lyapunov_candidate( + V_expr="x**2", + state_transition={"x": "2*x"}, + state_symbols=["x"], + continuous=False, + ) + assert result.stable is False + + def test_details_populated(self) -> None: + result = lyapunov_candidate( + V_expr="x**2", + state_transition={"x": "-x"}, + state_symbols=["x"], + ) + assert len(result.details) > 0 + assert result.dV_expr is not None + + +class TestQuadraticLyapunov: + def test_stable_diagonal(self) -> None: + """A = [[-1, 0], [0, -2]], P = I: A'P + PA = diag(-2, -4) < 0.""" + result = quadratic_lyapunov( + P=[[1.0, 0.0], [0.0, 1.0]], + A=[[-1.0, 0.0], [0.0, -2.0]], + ) + assert result.positive_definite == "PROVED" + assert result.decreasing == "PROVED" + assert result.stable is True + + def test_unstable_system(self) -> None: + """A = [[1, 0], [0, -1]], P = I: A'P + PA = diag(2, -2) not neg def.""" + result = quadratic_lyapunov( + P=[[1.0, 0.0], [0.0, 1.0]], + A=[[1.0, 0.0], [0.0, -1.0]], + ) + assert result.decreasing == "FAILED" + assert result.stable is False + + def test_custom_state_symbols(self) -> None: + result = quadratic_lyapunov( + P=[[1.0, 0.0], [0.0, 1.0]], + A=[[-1.0, 0.0], [0.0, -2.0]], + state_symbols=["pos", "vel"], + ) + assert result.stable is True + assert "pos" in str(result.candidate) + + +class TestFindQuadraticLyapunov: + def test_stable_system(self) -> None: + """Should find P for a stable system.""" + result = find_quadratic_lyapunov(A=[[-1.0, 0.0], [0.0, -2.0]]) + assert result is not None + P, lyap_result = result + assert lyap_result.stable is True + # P should be 2x2 + assert len(P) == 2 + assert len(P[0]) == 2 + + def test_unstable_system(self) -> None: + """Should return None for an unstable system.""" + result = find_quadratic_lyapunov(A=[[1.0, 0.0], [0.0, -1.0]]) + assert result is None + + def test_with_custom_Q(self) -> None: + result = find_quadratic_lyapunov( + A=[[-1.0, 0.0], [0.0, -2.0]], + Q=[[2.0, 0.0], [0.0, 2.0]], + ) + assert result is not None + + def test_empty_system(self) -> None: + assert find_quadratic_lyapunov(A=[]) is None + + +class TestPassivityCertificate: + def test_damped_spring(self) -> None: + """Spring-mass-damper: V = 0.5*(k*x^2 + m*v^2), supply = u*v. + + dx/dt = v, dv/dt = (-k*x - b*v + u) / m + With k=1, m=1, b=1: + dV/dt = k*x*v + m*v*(-k*x - b*v + u)/m + = k*x*v + v*(-k*x - b*v + u) + = k*x*v - k*x*v - b*v^2 + u*v + = -b*v^2 + u*v + Supply rate = u*v + Dissipation = dV/dt - supply = -b*v^2 <= 0 ✓ + """ + result = passivity_certificate( + V_expr="x**2/2 + v**2/2", + supply_rate="u*y_out", + state_transition={"x": "v", "v": "-x - v + u"}, + output_map={"y_out": "v"}, + state_symbols=["x", "v"], + input_symbols=["u"], + ) + + assert result.dissipation_proved == "PROVED" + assert result.passive is True + + def test_details_populated(self) -> None: + result = passivity_certificate( + V_expr="x**2", + supply_rate="u*y", + state_transition={"x": "-x + u"}, + output_map={"y": "x"}, + state_symbols=["x"], + input_symbols=["u"], + ) + assert len(result.details) > 0 diff --git a/packages/gds-viz/gds_viz/__init__.py b/packages/gds-viz/gds_viz/__init__.py index 4fda7e3..8a4147b 100644 --- a/packages/gds-viz/gds_viz/__init__.py +++ b/packages/gds-viz/gds_viz/__init__.py @@ -19,10 +19,28 @@ ] +_PHASE_EXPORTS = {"phase_portrait"} + +_FREQUENCY_EXPORTS = {"bode_plot", "nyquist_plot", "nichols_plot", "root_locus_plot"} + +_RESPONSE_EXPORTS = {"step_response_plot", "impulse_response_plot", "compare_responses"} + + def __getattr__(name: str) -> object: - """Lazy import for optional phase portrait module.""" - if name == "phase_portrait": + """Lazy import for optional visualization modules.""" + if name in _PHASE_EXPORTS: from gds_viz.phase import phase_portrait return phase_portrait + + if name in _FREQUENCY_EXPORTS: + from gds_viz import frequency + + return getattr(frequency, name) + + if name in _RESPONSE_EXPORTS: + from gds_viz import response + + return getattr(response, name) + raise AttributeError(f"module 'gds_viz' has no attribute {name!r}") diff --git a/packages/gds-viz/gds_viz/frequency.py b/packages/gds-viz/gds_viz/frequency.py new file mode 100644 index 0000000..3dc5daf --- /dev/null +++ b/packages/gds-viz/gds_viz/frequency.py @@ -0,0 +1,455 @@ +"""Frequency response visualization: Bode, Nyquist, Nichols, root locus. + +All functions accept plain ``list[float]`` or numpy arrays as input data. +No imports from gds-domains or gds-analysis — the caller computes the data, +the viz just plots it. + +Requires ``gds-viz[control]`` (matplotlib + numpy). +""" + +from __future__ import annotations + +from typing import Any + + +def _require_control_deps() -> None: + """Raise ImportError if matplotlib/numpy are absent.""" + try: + import matplotlib # noqa: F401 + import numpy # noqa: F401 + except ImportError as exc: + raise ImportError( + "Frequency response visualization requires matplotlib and numpy. " + "Install with: uv add gds-viz[control]" + ) from exc + + +def bode_plot( + omega: list[float], + mag_db: list[float], + phase_deg: list[float], + *, + title: str = "Bode Plot", + gain_margin: tuple[float, float] | None = None, + phase_margin: tuple[float, float] | None = None, + ax: tuple[Any, Any] | None = None, + figsize: tuple[float, float] = (10, 8), +) -> Any: + """Bode magnitude and phase plot. + + Two vertically stacked subplots: + - Top: magnitude (dB) vs. frequency (rad/s, log scale) + - Bottom: phase (degrees) vs. frequency (rad/s, log scale) + + Parameters + ---------- + omega : list[float] + Frequency values in rad/s. + mag_db : list[float] + Magnitude in dB. + phase_deg : list[float] + Phase in degrees. + title : str + Plot title. + gain_margin : tuple[float, float] | None + If provided, ``(gm_db, gm_freq)`` — annotates gain margin. + phase_margin : tuple[float, float] | None + If provided, ``(pm_deg, pm_freq)`` — annotates phase margin. + ax : tuple[Axes, Axes] | None + Existing (mag_ax, phase_ax) to plot into. Creates new if None. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + w = np.array(omega) + mag = np.array(mag_db) + phase = np.array(phase_deg) + + if ax is None: + fig, (ax_mag, ax_phase) = plt.subplots(2, 1, figsize=figsize, sharex=True) + else: + ax_mag, ax_phase = ax + fig = ax_mag.get_figure() + + # Magnitude plot + ax_mag.semilogx(w, mag, "b-", linewidth=1.5) + ax_mag.axhline(y=0, color="k", linewidth=0.5, linestyle="--", alpha=0.5) + ax_mag.set_ylabel("Magnitude (dB)") + ax_mag.set_title(title) + ax_mag.grid(True, which="both", alpha=0.3) + + # Phase plot + ax_phase.semilogx(w, phase, "r-", linewidth=1.5) + ax_phase.axhline(y=-180, color="k", linewidth=0.5, linestyle="--", alpha=0.5) + ax_phase.set_xlabel("Frequency (rad/s)") + ax_phase.set_ylabel("Phase (deg)") + ax_phase.grid(True, which="both", alpha=0.3) + + # Annotate margins + if gain_margin is not None: + gm_db, gm_freq = gain_margin + if np.isfinite(gm_db) and np.isfinite(gm_freq): + ax_mag.axvline(x=gm_freq, color="g", linestyle=":", alpha=0.7) + ax_mag.annotate( + f"GM = {gm_db:.1f} dB", + xy=(gm_freq, 0), + xytext=(gm_freq * 2, 5), + fontsize=9, + color="g", + arrowprops={"arrowstyle": "->", "color": "g"}, + ) + + if phase_margin is not None: + pm_deg_val, pm_freq = phase_margin + if np.isfinite(pm_deg_val) and np.isfinite(pm_freq): + ax_phase.axvline(x=pm_freq, color="m", linestyle=":", alpha=0.7) + ax_phase.annotate( + f"PM = {pm_deg_val:.1f}°", + xy=(pm_freq, -180), + xytext=(pm_freq * 2, -180 + pm_deg_val / 2), + fontsize=9, + color="m", + arrowprops={"arrowstyle": "->", "color": "m"}, + ) + + plt.tight_layout() + return fig + + +def nyquist_plot( + real: list[float], + imag: list[float], + *, + title: str = "Nyquist Plot", + unit_circle: bool = True, + critical_point: bool = True, + ax: Any | None = None, + figsize: tuple[float, float] = (8, 8), +) -> Any: + """Nyquist diagram (polar frequency response). + + Plots H(jω) in the complex plane. + + Parameters + ---------- + real : list[float] + Real part of H(jω). + imag : list[float] + Imaginary part of H(jω). + title : str + Plot title. + unit_circle : bool + If True, draw the unit circle. + critical_point : bool + If True, mark the critical point (-1, 0). + ax : Axes | None + Existing axes. Creates new if None. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + re = np.array(real) + im = np.array(imag) + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize) + else: + fig = ax.get_figure() + + # Positive frequencies + ax.plot(re, im, "b-", linewidth=1.5, label="ω > 0") + # Mirror (negative frequencies) + ax.plot(re, -im, "b--", linewidth=0.8, alpha=0.5, label="ω < 0") + + if unit_circle: + theta = np.linspace(0, 2 * np.pi, 200) + ax.plot(np.cos(theta), np.sin(theta), "k--", linewidth=0.5, alpha=0.3) + + if critical_point: + ax.plot(-1, 0, "rx", markersize=10, markeredgewidth=2, label="(-1, 0)") + + # Mark start and end + ax.plot(re[0], im[0], "go", markersize=6, label=f"ω={0}") + if len(re) > 1: + ax.annotate( + "", + xy=(re[len(re) // 2], im[len(re) // 2]), + xytext=(re[len(re) // 2 - 1], im[len(re) // 2 - 1]), + arrowprops={"arrowstyle": "->", "color": "b"}, + ) + + ax.set_xlabel("Real") + ax.set_ylabel("Imaginary") + ax.set_title(title) + ax.set_aspect("equal") + ax.grid(True, alpha=0.3) + ax.legend(fontsize=8) + plt.tight_layout() + return fig + + +def nichols_plot( + phase_deg: list[float], + mag_db: list[float], + *, + title: str = "Nichols Chart", + m_circles: bool = True, + n_circles: bool = False, + ax: Any | None = None, + figsize: tuple[float, float] = (10, 8), +) -> Any: + """Nichols chart: open-loop phase (x) vs. open-loop gain (y). + + Parameters + ---------- + phase_deg : list[float] + Open-loop phase in degrees. + mag_db : list[float] + Open-loop magnitude in dB. + title : str + Plot title. + m_circles : bool + If True, draw M-circles (constant closed-loop magnitude contours). + n_circles : bool + If True, draw N-circles (constant closed-loop phase contours). + ax : Axes | None + Existing axes. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + phase = np.array(phase_deg) + mag = np.array(mag_db) + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize) + else: + fig = ax.get_figure() + + ax.plot(phase, mag, "b-", linewidth=1.5) + ax.plot(phase[0], mag[0], "go", markersize=6) + ax.plot(phase[-1], mag[-1], "rs", markersize=5) + + # Mark critical point (-180°, 0 dB) + ax.plot(-180, 0, "rx", markersize=10, markeredgewidth=2) + + if m_circles: + # Draw M-circles at standard dB levels + m_values_db = [-12, -6, -3, -1, 0.5, 1, 3, 6, 12] + theta = np.linspace(-359, 0, 1000) + + for m_db in m_values_db: + M = 10 ** (m_db / 20) + if abs(M - 1.0) < 1e-10: + continue + # M-circle in Nichols coordinates + # |T| = M where T = L/(1+L), L = G*e^{j*phase} + # Parametric: for each open-loop phase, compute the gain + # that gives |T| = M + phase_rad = np.radians(theta) + # |L / (1+L)| = M => |L|^2 = M^2 |1+L|^2 + # Let g = |L|: g^2 = M^2(1 + 2g*cos(phi) + g^2) + # g^2(1-M^2) = M^2(1 + 2g*cos(phi)) + # This is quadratic in g + cos_phi = np.cos(phase_rad) + a_coeff = 1 - M**2 + b_coeff = -2 * M**2 * cos_phi + c_coeff = -(M**2) + + discriminant = b_coeff**2 - 4 * a_coeff * c_coeff + valid = discriminant >= 0 + g = np.full_like(theta, np.nan) + if abs(a_coeff) > 1e-10: + g[valid] = (-b_coeff[valid] + np.sqrt(discriminant[valid])) / ( + 2 * a_coeff + ) + g_db = 20 * np.log10(np.maximum(g, 1e-15)) + mask = (g > 0) & np.isfinite(g_db) & (np.abs(g_db) < 40) + if np.any(mask): + ax.plot( + theta[mask], + g_db[mask], + "k-", + linewidth=0.3, + alpha=0.3, + ) + + ax.set_xlabel("Open-Loop Phase (deg)") + ax.set_ylabel("Open-Loop Gain (dB)") + ax.set_title(title) + ax.set_xlim(-360, 0) + ax.grid(True, alpha=0.3) + plt.tight_layout() + return fig + + +def root_locus_plot( + num: list[float], + den: list[float], + *, + gains: list[float] | None = None, + title: str = "Root Locus", + mark_poles: bool = True, + mark_zeros: bool = True, + ax: Any | None = None, + figsize: tuple[float, float] = (8, 8), +) -> Any: + """Root locus diagram. + + Shows pole migration as gain K varies from 0 to max. + + Parameters + ---------- + num : list[float] + Transfer function numerator coefficients (descending powers). + den : list[float] + Transfer function denominator coefficients (descending powers). + gains : list[float] | None + Gain values to evaluate. If None, auto-selects 500 values. + title : str + Plot title. + mark_poles : bool + If True, mark open-loop poles with x. + mark_zeros : bool + If True, mark open-loop zeros with o. + ax : Axes | None + Existing axes. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + num_arr = np.array(num, dtype=float) + den_arr = np.array(den, dtype=float) + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize) + else: + fig = ax.get_figure() + + if gains is None: + # Auto-select gains: logarithmic sweep + gains_arr = np.concatenate( + [ + np.array([0.0]), + np.logspace(-3, 3, 500), + ] + ) + else: + gains_arr = np.array(gains, dtype=float) + + # Compute closed-loop poles for each gain + n_poles = len(den_arr) - 1 + all_poles = np.zeros((len(gains_arr), n_poles), dtype=complex) + + for i, k in enumerate(gains_arr): + # Closed-loop characteristic polynomial: den + K * num + # Pad num to same length as den + num_padded = np.zeros(len(den_arr)) + num_padded[-len(num_arr) :] = num_arr + cl_poly = den_arr + k * num_padded + roots = np.roots(cl_poly) + + # Sort roots to maintain continuity + if i > 0: + roots = _sort_poles_by_proximity(all_poles[i - 1], roots) + all_poles[i] = roots + + # Plot pole trajectories + cmap = plt.get_cmap("viridis") + for j in range(n_poles): + trajectory = all_poles[:, j] + colors = [cmap(i / max(len(gains_arr) - 1, 1)) for i in range(len(gains_arr))] + for k in range(1, len(trajectory)): + ax.plot( + [trajectory[k - 1].real, trajectory[k].real], + [trajectory[k - 1].imag, trajectory[k].imag], + "-", + color=colors[k], + linewidth=0.8, + alpha=0.7, + ) + + # Mark open-loop poles and zeros + if mark_poles: + ol_poles = np.roots(den_arr) + ax.plot( + ol_poles.real, + ol_poles.imag, + "x", + color="red", + markersize=10, + markeredgewidth=2, + label="Open-loop poles", + ) + + if mark_zeros and len(num_arr) > 1: + ol_zeros = np.roots(num_arr) + ax.plot( + ol_zeros.real, + ol_zeros.imag, + "o", + color="blue", + markersize=8, + markerfacecolor="none", + markeredgewidth=2, + label="Open-loop zeros", + ) + + # Imaginary axis + ax.axvline(x=0, color="k", linewidth=0.5, alpha=0.3) + ax.axhline(y=0, color="k", linewidth=0.5, alpha=0.3) + + ax.set_xlabel("Real") + ax.set_ylabel("Imaginary") + ax.set_title(title) + ax.grid(True, alpha=0.3) + ax.legend(fontsize=8) + ax.set_aspect("equal") + plt.tight_layout() + return fig + + +def _sort_poles_by_proximity(prev_poles: Any, new_poles: Any) -> Any: + """Sort new_poles to minimize distance from prev_poles (greedy matching).""" + import numpy as np + + result = np.zeros_like(new_poles) + available = list(range(len(new_poles))) + + for i in range(len(prev_poles)): + if not available: + break + distances = [abs(new_poles[j] - prev_poles[i]) for j in available] + best_idx = available[int(np.argmin(distances))] + result[i] = new_poles[best_idx] + available.remove(best_idx) + + return result diff --git a/packages/gds-viz/gds_viz/response.py b/packages/gds-viz/gds_viz/response.py new file mode 100644 index 0000000..fcac33c --- /dev/null +++ b/packages/gds-viz/gds_viz/response.py @@ -0,0 +1,272 @@ +"""Step and impulse response visualization. + +Provides annotated step response plots with performance metric annotations +(rise time, settling time, overshoot bands) and multi-response comparison. + +All functions accept plain ``list[float]`` data. No imports from +gds-domains or gds-analysis — the caller computes data, the viz plots it. + +Requires ``gds-viz[control]`` (matplotlib + numpy). +""" + +from __future__ import annotations + +from typing import Any + + +def _require_control_deps() -> None: + """Raise ImportError if matplotlib/numpy are absent.""" + try: + import matplotlib # noqa: F401 + import numpy # noqa: F401 + except ImportError as exc: + raise ImportError( + "Response visualization requires matplotlib and numpy. " + "Install with: uv add gds-viz[control]" + ) from exc + + +def step_response_plot( + times: list[float], + values: list[float], + *, + setpoint: float | None = None, + metrics: Any | None = None, + title: str = "Step Response", + ax: Any | None = None, + figsize: tuple[float, float] = (10, 6), +) -> Any: + """Step response plot with optional metric annotations. + + Parameters + ---------- + times : list[float] + Time values. + values : list[float] + Response values. + setpoint : float | None + Desired final value. If provided, draws a reference line. + metrics : StepMetrics | None + If provided (from ``gds_analysis.response``), annotates: + - Rise time (shaded region) + - Settling time (vertical dashed line) + - Overshoot (horizontal dashed line at peak) + - ±2% settling band (shaded horizontal band) + title : str + Plot title. + ax : Axes | None + Existing axes. Creates new if None. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + t = np.array(times) + y = np.array(values) + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize) + else: + fig = ax.get_figure() + + ax.plot(t, y, "b-", linewidth=1.5, label="Response") + + if setpoint is not None: + ax.axhline( + y=setpoint, + color="k", + linewidth=0.8, + linestyle="--", + alpha=0.5, + label=f"Setpoint = {setpoint}", + ) + + if metrics is not None: + _annotate_metrics(ax, t, y, metrics) + + ax.set_xlabel("Time (s)") + ax.set_ylabel("Response") + ax.set_title(title) + ax.grid(True, alpha=0.3) + ax.legend(fontsize=8) + plt.tight_layout() + return fig + + +def _annotate_metrics(ax: Any, t: Any, y: Any, metrics: Any) -> None: + """Add metric annotations to a step response plot.""" + ss = metrics.steady_state_value + + # Settling band (±2%) + band = 0.02 * abs(ss) if abs(ss) > 1e-15 else 0.02 + ax.axhspan( + ss - band, + ss + band, + color="green", + alpha=0.08, + label="±2% band", + ) + + # Settling time + if metrics.settling_time > t[0]: + ax.axvline( + x=metrics.settling_time, + color="green", + linewidth=1, + linestyle="--", + alpha=0.7, + label=f"Settling: {metrics.settling_time:.2f}s", + ) + + # Rise time + if metrics.rise_time > 0: + v0 = y[0] + v_range = ss - v0 + t_low = _find_crossing_time(t, y, v0 + 0.1 * v_range) + t_high = _find_crossing_time(t, y, v0 + 0.9 * v_range) + if t_low < t_high: + ax.axvspan( + t_low, + t_high, + color="orange", + alpha=0.1, + label=f"Rise: {metrics.rise_time:.2f}s", + ) + + # Overshoot + if metrics.overshoot_pct > 0: + ax.axhline( + y=metrics.peak_value, + color="red", + linewidth=0.8, + linestyle=":", + alpha=0.6, + ) + ax.plot( + metrics.peak_time, + metrics.peak_value, + "rv", + markersize=8, + label=f"Overshoot: {metrics.overshoot_pct:.1f}%", + ) + + # Steady-state error + if metrics.steady_state_error > 0.001: + ax.annotate( + f"SS error: {metrics.steady_state_error:.4f}", + xy=(t[-1] * 0.8, ss), + fontsize=8, + color="gray", + ) + + +def _find_crossing_time(t: Any, y: Any, threshold: float) -> float: + """Find first time where y crosses threshold.""" + for i in range(1, len(y)): + if (y[i - 1] <= threshold <= y[i]) or (y[i - 1] >= threshold >= y[i]): + dv = y[i] - y[i - 1] + if abs(dv) < 1e-15: + return float(t[i]) + frac = (threshold - y[i - 1]) / dv + return float(t[i - 1] + frac * (t[i] - t[i - 1])) + return float(t[-1]) + + +def impulse_response_plot( + times: list[float], + values: list[float], + *, + title: str = "Impulse Response", + ax: Any | None = None, + figsize: tuple[float, float] = (10, 6), +) -> Any: + """Impulse response plot. + + Parameters + ---------- + times : list[float] + Time values. + values : list[float] + Impulse response values. + title : str + Plot title. + ax : Axes | None + Existing axes. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + t = np.array(times) + y = np.array(values) + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize) + else: + fig = ax.get_figure() + + ax.plot(t, y, "b-", linewidth=1.5) + ax.axhline(y=0, color="k", linewidth=0.5, alpha=0.3) + ax.set_xlabel("Time (s)") + ax.set_ylabel("Response") + ax.set_title(title) + ax.grid(True, alpha=0.3) + plt.tight_layout() + return fig + + +def compare_responses( + responses: list[tuple[list[float], list[float], str]], + *, + title: str = "Response Comparison", + figsize: tuple[float, float] = (10, 6), +) -> Any: + """Overlay multiple step/impulse responses for comparison. + + Parameters + ---------- + responses : list[tuple[list[float], list[float], str]] + Each tuple is ``(times, values, label)``. + title : str + Plot title. + figsize : tuple[float, float] + Figure size. + + Returns + ------- + matplotlib Figure + """ + _require_control_deps() + import matplotlib.pyplot as plt + import numpy as np + + fig, ax = plt.subplots(1, 1, figsize=figsize) + + cmap = plt.get_cmap("tab10") + n = max(len(responses), 1) + + for i, (times, values, label) in enumerate(responses): + t = np.array(times) + y = np.array(values) + color = cmap(i / n) + ax.plot(t, y, "-", color=color, linewidth=1.5, label=label) + + ax.set_xlabel("Time (s)") + ax.set_ylabel("Response") + ax.set_title(title) + ax.grid(True, alpha=0.3) + ax.legend(fontsize=8) + plt.tight_layout() + return fig diff --git a/packages/gds-viz/pyproject.toml b/packages/gds-viz/pyproject.toml index 7b718d1..7e6f56a 100644 --- a/packages/gds-viz/pyproject.toml +++ b/packages/gds-viz/pyproject.toml @@ -35,6 +35,8 @@ dependencies = [ [project.optional-dependencies] phase = ["matplotlib>=3.8", "numpy>=1.26", "gds-continuous>=0.1.0"] +control = ["matplotlib>=3.8", "numpy>=1.26"] +plots = ["matplotlib>=3.8", "numpy>=1.26", "gds-continuous>=0.1.0"] [project.urls] Homepage = "https://github.com/BlockScience/gds-core" From 0fc24dbf5fb6d1087dd1d02f43f823c9d3618ad7 Mon Sep 17 00:00:00 2001 From: rohan Date: Tue, 7 Apr 2026 19:01:16 +0530 Subject: [PATCH 2/3] fix: test coverage, docs, and cleanup for control theory stack MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Quick wins from the post-implementation audit: Tests: - Add 26 viz tests: Bode, Nyquist, Nichols, root locus, step/impulse response plots, compare_responses, metric annotations - Add gain_margin numerical verification (35.6 dB for third-order system) - Add phase_margin with finite margin (integrator: PM = 90 deg) - Add MIMO transfer function tests (2x2, DC gains, nonzero D feedthrough) - Add Gang of Six DC value verification (S(0), T(0), CS(0)) - Add gain_schedule test (spring-mass at 3 operating points) - Add LQR with cross-term N test - Add backward_euler discretization test Fixes: - Remove unimplemented n_circles parameter from nichols_plot - Fix stale docstring in lyapunov.py (no longer claims to delegate to the five-strategy prover) - Document _find_crossing_time duplication between gds-viz and gds-analysis (intentional: avoids hard dependency) Docs: - Add CHANGELOG entry for the full control theory stack - Update root CLAUDE.md with symbolic→analysis→viz data flow - Update gds-symbolic CLAUDE.md with transfer.py and delay.py - Update gds-analysis CLAUDE.md with linear.py and response.py - Update gds-viz CLAUDE.md with frequency/response plot sections - Update gds-proof CLAUDE.md with Lyapunov module documentation --- CHANGELOG.md | 65 +++++++ CLAUDE.md | 20 ++- packages/gds-analysis/CLAUDE.md | 20 ++- packages/gds-analysis/tests/test_linear.py | 160 +++++++++++++++--- packages/gds-domains/tests/test_transfer.py | 119 ++++++++++++- packages/gds-proof/CLAUDE.md | 18 ++ .../gds-proof/gds_proof/analysis/lyapunov.py | 4 +- packages/gds-symbolic/CLAUDE.md | 17 ++ packages/gds-viz/CLAUDE.md | 21 ++- packages/gds-viz/gds_viz/frequency.py | 3 - packages/gds-viz/gds_viz/response.py | 6 +- packages/gds-viz/tests/test_frequency.py | 129 ++++++++++++++ packages/gds-viz/tests/test_response_viz.py | 140 +++++++++++++++ 13 files changed, 677 insertions(+), 45 deletions(-) create mode 100644 packages/gds-viz/tests/test_frequency.py create mode 100644 packages/gds-viz/tests/test_response_viz.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 271ab8a..c308756 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,70 @@ # Changelog +## 2026-04-07 — Classical Control Theory Analysis Stack + +Adds frequency-domain and time-domain control analysis across four packages, +closing the gap between GDS structural specification and classical control +engineering workflows. Tracking issue: #198. + +### gds-domains (symbolic layer) + +- **`symbolic/transfer.py`** — transfer function representation and analysis + - `ss_to_tf()`: state-space to transfer function matrix (SISO and MIMO) + - `poles()`, `zeros()`: polynomial root-finding via SymPy + - `characteristic_polynomial()`: det(sI - A) as coefficient list + - `controllability_matrix()`, `observability_matrix()`: Kalman rank tests + - `is_controllable()`, `is_observable()`: rank-based checks + - `is_minimum_phase()`: RHP zero detection + - `sensitivity()`: Gang of Six (S, T, CS, PS, KS, KPS) + - `TransferFunction`, `TransferFunctionMatrix` data types +- **`symbolic/delay.py`** — time delay modeling + - `pade_approximation()`: (N,N) Padé approximant of e^{-sτ} + - `delay_system()`: cascade TF with Padé delay + +### gds-analysis (numerical layer) + +- **`linear.py`** — numerical linear systems analysis (requires `[continuous]`) + - `eigenvalues()`, `is_stable()`, `is_marginally_stable()`: stability checks + - `frequency_response()`: H(jω) magnitude/phase via scipy + - `gain_margin()`, `phase_margin()`: stability margin computation + - `discretize()`: continuous-to-discrete (Tustin, ZOH, Euler, backward Euler) + - `lqr()`, `dlqr()`: continuous/discrete Linear Quadratic Regulator + - `kalman()`: steady-state Kalman filter gain + - `gain_schedule()`: multi-point LQR across operating points +- **`response.py`** — step/impulse response and time-domain metrics + - `step_response()`, `impulse_response()`: state-space simulation via scipy + - `step_response_metrics()`: rise time, settling time, overshoot, SSE + - `StepMetrics` data type (no scipy needed for metric extraction) + - `metrics_from_ode_results()`: convenience wrapper for gds-continuous + +### gds-proof (formal verification layer) + +- **`analysis/lyapunov.py`** — Lyapunov stability proofs + - `lyapunov_candidate()`: verify V(x) > 0 and dV/dt < 0 (continuous/discrete) + - `quadratic_lyapunov()`: verify V = x'Px via A'P + PA eigenvalue check + - `find_quadratic_lyapunov()`: solve Lyapunov equation A'P + PA = -Q + - `passivity_certificate()`: verify dV/dt ≤ s(u, y) for dissipativity + +### gds-viz (visualization layer) + +- **`frequency.py`** — frequency response plots (requires `[control]`) + - `bode_plot()`: magnitude + phase with optional margin annotations + - `nyquist_plot()`: complex plane with unit circle and critical point + - `nichols_plot()`: open-loop phase vs. gain with M-circles + - `root_locus_plot()`: pole migration with gain sweep +- **`response.py`** — time-domain plots (requires `[control]`) + - `step_response_plot()`: with StepMetrics annotation overlays + - `impulse_response_plot()`: basic response plot + - `compare_responses()`: multi-response overlay +- New `[control]` optional extra (matplotlib + numpy, no gds-continuous) +- New `[plots]` extra (union of `[phase]` and `[control]`) + +Zero new third-party dependencies. 104 new tests across all packages. + +Closes #199, #200, #201, #202, #203, #204. + +--- + ## 2026-04-07 — gds-proof Layer 1 Integration Promotes gds-proof from a standalone package (Layer 0, protocol-only) to a diff --git a/CLAUDE.md b/CLAUDE.md index 444e500..04d7969 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -14,13 +14,27 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co | gds-domains | `gds_domains.{stockflow,control,business,software,games,symbolic}` | `packages/gds-domains/` | All domain DSLs | | gds-sim | `gds_sim` | `packages/gds-sim/` | Discrete-time simulation | | gds-continuous | `gds_continuous` | `packages/gds-continuous/` | Continuous-time ODE engine | -| gds-analysis | `gds_analysis`, `gds_analysis.psuu` | `packages/gds-analysis/` | Spec-to-sim bridge + PSUU | +| gds-analysis | `gds_analysis`, `gds_analysis.psuu` | `packages/gds-analysis/` | Spec-to-sim bridge + PSUU + linear systems analysis | | gds-interchange | `gds_interchange.owl` | `packages/gds-interchange/` | OWL/SHACL/SPARQL + future bridges | -| gds-viz | `gds_viz` | `packages/gds-viz/` | Mermaid + phase portraits | +| gds-viz | `gds_viz` | `packages/gds-viz/` | Mermaid + phase portraits + Bode/Nyquist/root locus | | gds-examples | — | `packages/gds-examples/` | Tutorials + examples | Deprecated shim packages (v0.99.0, re-export with DeprecationWarning): gds-owl, gds-psuu, gds-stockflow, gds-control, gds-games, gds-software, gds-business, gds-symbolic. +## What GDS Is Not + +GDS is a **specification and verification** framework, not a simulation engine or general-purpose modeling library. If you have used SimPy, Mesa, cadCAD, or SciPy before, read this section before writing any code. + +| Misconception | Reality | +|--------------|---------| +| GDS is a simulation engine | GDS is a **specification** layer. Do not write `for t in range(...)` loops. Declare a spec; `gds-sim` or `gds-continuous` executes it. | +| GDS is SimPy / Mesa / cadCAD | GDS is a typed composition algebra with formal verification. The closest analogy is a specification language (like Alloy or TLA+), not an ABM or process-simulation framework. | +| GDS is SciPy.integrate | `gds-continuous` wraps SciPy internally, but users declare ODE systems via specs, not `solve_ivp()` calls directly. | +| Blocks are classes you instantiate freely | Block roles (`BoundaryAction`, `Policy`, `Mechanism`, `ControlAction`) are the four leaf types. Domain DSLs compile *models* into blocks via `compile_model()` / `compile_to_system()`. Hand-wiring blocks is for framework extension, not typical usage. | +| Verification is type checking | Verification operates on compiled IR (`SystemIR`) and specs (`GDSSpec`). It checks structural topology (G-001..G-006) and semantic properties (SC-001..SC-009), not Python types at runtime. | +| Parameters get assigned values by GDS | `ParameterSchema` (Theta) is structural metadata only. GDS never binds parameter values. `gds-analysis` and its PSUU module handle sweep and evaluation. | +| Domain DSLs are interchangeable | Each DSL targets a specific modeling paradigm. Use `gds_domains.stockflow` for conservation-law systems, `.control` for state-space controllers, `.games` for strategic interaction, `.business` for causal/value-stream diagrams, `.software` for architecture diagrams. See the [Choosing a DSL](https://blockscience.github.io/gds-core/guides/choosing-a-dsl/) guide. | + ## Commands ```bash @@ -107,6 +121,8 @@ The composition tree follows a convergent tiered pattern: `gds_domains.games` is more complex — it subclasses `AtomicBlock` as `OpenGame` with its own IR (`PatternIR`), but projects back to `SystemIR` via `PatternIR.to_system_ir()`. +`gds_domains.symbolic` extends the control DSL with SymPy-based ODEs, linearization (A, B, C, D matrices), transfer functions, poles/zeros, controllability/observability, Padé time delay approximation, and the Gang of Six sensitivity functions. This is the entry point to the classical control theory analysis stack — `LinearizedSystem` flows into `gds-analysis` for numerical analysis and `gds-viz` for frequency-domain plots. + ### Two Type Systems These coexist at different levels: diff --git a/packages/gds-analysis/CLAUDE.md b/packages/gds-analysis/CLAUDE.md index f999ac1..88fb3d0 100644 --- a/packages/gds-analysis/CLAUDE.md +++ b/packages/gds-analysis/CLAUDE.md @@ -10,7 +10,7 @@ ## Architecture -Five modules, each bridging one aspect of structural specification to runtime: +Seven modules bridging structural specification to runtime analysis: | Module | Function | Paper reference | |--------|----------|-----------------| @@ -19,6 +19,8 @@ Five modules, each bridging one aspect of structural specification to runtime: | `metrics.py` | `trajectory_distances(results, spec)` → distance matrix | Assumption 3.2 | | `reachability.py` | `reachable_set(model, state, inputs)` → R(x) | Def 4.1, 4.2 | | `backward_reachability.py` | `backward_reachable_set(dynamics, ...)` → B(T) | Def 4.1 (backward) | +| `linear.py` | Eigenvalue stability, frequency response, margins, discretization, LQR, Kalman | `[continuous]` | +| `response.py` | Step/impulse response computation + time-domain metrics (StepMetrics) | `[continuous]` | ### spec_to_model adapter @@ -35,6 +37,22 @@ If `enforce_constraints=True`, wraps BoundaryAction policies with `guarded_polic - `reachable_graph(spec, model, states, input_samples)` — builds full reachability graph across multiple states - `configuration_space(reachability_graph)` — extracts largest SCC (the configuration space X_C) +### Linear systems analysis (`linear.py`, requires `[continuous]`) + +All functions accept `list[list[float]]` matrices (matching `LinearizedSystem` fields): +- `eigenvalues(A)`, `is_stable(A)`, `is_marginally_stable(A)` — stability checks +- `frequency_response(A, B, C, D, omega)` → `(omega, mag_db, phase_deg)` +- `gain_margin(num, den)`, `phase_margin(num, den)` — stability margins +- `discretize(A, B, C, D, dt, method)` → `(Ad, Bd, Cd, Dd)` via scipy.signal +- `lqr(A, B, Q, R)`, `dlqr(Ad, Bd, Q, R)` → `(K, P, E)` gain + Riccati + eigenvalues +- `kalman(A, C, Q, R)` → `(L, P)` observer gain + covariance +- `gain_schedule(linearize_fn, points, Q, R)` → gains at multiple operating points + +### Response metrics (`response.py`) + +- `step_response_metrics(times, values, setpoint)` → `StepMetrics` (no scipy needed) +- `step_response(A, B, C, D)`, `impulse_response(A, B, C, D)` — scipy-based simulation + ## Commands ```bash diff --git a/packages/gds-analysis/tests/test_linear.py b/packages/gds-analysis/tests/test_linear.py index 084ec34..78ed64d 100644 --- a/packages/gds-analysis/tests/test_linear.py +++ b/packages/gds-analysis/tests/test_linear.py @@ -8,6 +8,7 @@ eigenvalues, frequency_response, gain_margin, + gain_schedule, is_marginally_stable, is_stable, kalman, @@ -15,7 +16,6 @@ phase_margin, ) - # --------------------------------------------------------------------------- # Stability analysis # --------------------------------------------------------------------------- @@ -85,7 +85,10 @@ class TestFrequencyResponse: def test_first_order_bode(self) -> None: """1/(s+1): -3dB at w=1, -45deg at w=1.""" omega, mag_db, phase_deg = frequency_response( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], omega=[1.0], ) @@ -94,8 +97,11 @@ def test_first_order_bode(self) -> None: assert phase_deg[0] == pytest.approx(-45.0, abs=2.0) def test_auto_frequency_range(self) -> None: - omega, mag_db, phase_deg = frequency_response( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], + omega, _mag_db, _phase_deg = frequency_response( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], ) assert len(omega) == 500 # default n_points @@ -106,20 +112,43 @@ def test_known_system(self) -> None: # H(s) = 1/(s^3 + 6s^2 + 11s + 6) num = [1.0] den = [1.0, 6.0, 11.0, 6.0] - gm_db, gm_freq = gain_margin(num, den) + gm_db, _gm_freq = gain_margin(num, den) assert gm_db > 0 # System is stable, should have positive margin + def test_third_order_numerical_value(self) -> None: + """1/((s+1)(s+2)(s+3)): verify gain margin and crossover frequency. + + Phase crosses -180 deg at w ≈ 3.32 rad/s. + |H(j*3.32)| = 1/sqrt((1+11)(4+11)(9+11)) ≈ 0.0167 → GM ≈ 35.6 dB. + """ + num = [1.0] + den = [1.0, 6.0, 11.0, 6.0] + gm_db, gm_freq = gain_margin(num, den) + assert gm_db == pytest.approx(35.6, abs=2.0) + assert gm_freq == pytest.approx(3.32, abs=0.3) + + def test_infinite_margin_stable_first_order(self) -> None: + """1/(s+1) never reaches -180° → infinite gain margin.""" + gm_db, _gm_freq = gain_margin([1.0], [1.0, 1.0]) + assert gm_db == float("inf") + class TestPhaseMargin: - def test_known_system(self) -> None: - """1/(s+1): infinite phase margin (gain never reaches 0dB at - a frequency where phase is problematic).""" + def test_infinite_for_first_order(self) -> None: + """1/(s+1): gain never reaches 0 dB → infinite phase margin.""" num = [1.0] den = [1.0, 1.0] - pm_deg, pm_freq = phase_margin(num, den) - # For 1/(s+1), gain is always < 0dB, so phase margin is infinite + pm_deg, _pm_freq = phase_margin(num, den) assert pm_deg == float("inf") + def test_integrator_with_gain(self) -> None: + """K/s with K=1: gain crosses 0 dB at w=1, phase = -90° → PM = 90°.""" + num = [1.0] + den = [1.0, 0.0] # 1/s + pm_deg, pm_freq = phase_margin(num, den) + assert pm_deg == pytest.approx(90.0, abs=5.0) + assert pm_freq == pytest.approx(1.0, abs=0.3) + # --------------------------------------------------------------------------- # Discretization @@ -129,35 +158,63 @@ def test_known_system(self) -> None: class TestDiscretize: def test_tustin_preserves_stability(self) -> None: """Stable continuous system → stable discrete system.""" - Ad, Bd, Cd, Dd = discretize( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], - dt=0.1, method="tustin", + Ad, _Bd, _Cd, _Dd = discretize( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + dt=0.1, + method="tustin", ) assert is_stable(Ad, continuous=False) def test_zoh(self) -> None: """ZOH discretization of 1/(s+1) with dt=0.1.""" - Ad, Bd, Cd, Dd = discretize( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], - dt=0.1, method="zoh", + Ad, _Bd, _Cd, _Dd = discretize( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + dt=0.1, + method="zoh", ) # A_d = e^{A*dt} = e^{-0.1} ≈ 0.9048 assert Ad[0][0] == pytest.approx(0.9048, abs=0.01) assert is_stable(Ad, continuous=False) def test_euler(self) -> None: - Ad, Bd, Cd, Dd = discretize( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], - dt=0.1, method="euler", + Ad, _Bd, _Cd, _Dd = discretize( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + dt=0.1, + method="euler", ) # Forward Euler: Ad = I + dt*A = [[0.9]] assert Ad[0][0] == pytest.approx(0.9, abs=0.01) + def test_backward_euler(self) -> None: + """Backward Euler preserves stability for any step size.""" + Ad, _Bd, _Cd, _Dd = discretize( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + dt=0.5, + method="backward_euler", + ) + assert is_stable(Ad, continuous=False) + def test_invalid_method(self) -> None: with pytest.raises(ValueError, match="Unknown discretization method"): discretize( - A=[[-1.0]], B=[[1.0]], C=[[1.0]], D=[[0.0]], - dt=0.1, method="invalid", + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[0.0]], + dt=0.1, + method="invalid", ) @@ -174,7 +231,7 @@ def test_double_integrator(self) -> None: Q = [[1.0, 0.0], [0.0, 1.0]] R = [[1.0]] - K, P, E = lqr(A, B, Q, R) + K, _P, E = lqr(A, B, Q, R) # K should be 1x2 assert len(K) == 1 @@ -185,30 +242,46 @@ def test_double_integrator(self) -> None: def test_scalar_system(self) -> None: """Simple scalar system: dx/dt = -x + u, Q=1, R=1.""" - K, P, E = lqr( - A=[[-1.0]], B=[[1.0]], Q=[[1.0]], R=[[1.0]], + K, _P, E = lqr( + A=[[-1.0]], + B=[[1.0]], + Q=[[1.0]], + R=[[1.0]], ) assert len(K) == 1 assert len(K[0]) == 1 assert E[0].real < 0 + def test_with_cross_term_N(self) -> None: + """LQR with cross-term N should still produce stable closed-loop.""" + A = [[0.0, 1.0], [0.0, 0.0]] + B = [[0.0], [1.0]] + Q = [[10.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + N = [[0.1], [0.0]] + + _K, _P, E = lqr(A, B, Q, R, N=N) + + assert all(e.real < 0 for e in E) + class TestDLQR: def test_discrete_double_integrator(self) -> None: """Discrete LQR on discretized double integrator.""" # First discretize - Ad, Bd, Cd, Dd = discretize( + Ad, Bd, _Cd, _Dd = discretize( A=[[0.0, 1.0], [0.0, 0.0]], B=[[0.0], [1.0]], C=[[1.0, 0.0]], D=[[0.0]], - dt=0.1, method="zoh", + dt=0.1, + method="zoh", ) Q = [[1.0, 0.0], [0.0, 1.0]] R = [[1.0]] - K, P, E = dlqr(Ad, Bd, Q, R) + _K, _P, E = dlqr(Ad, Bd, Q, R) # Closed-loop eigenvalues inside unit circle assert all(abs(e) < 1.0 for e in E) @@ -222,7 +295,7 @@ def test_observer_gain(self) -> None: Q_proc = [[1.0, 0.0], [0.0, 1.0]] R_meas = [[0.1]] - L, P = kalman(A, C, Q_proc, R_meas) + L, _P = kalman(A, C, Q_proc, R_meas) # L should be 2x1 (n x p) assert len(L) == 2 @@ -234,3 +307,34 @@ def test_observer_gain(self) -> None: A_obs = np.array(A) - np.array(L) @ np.array(C) eigs = np.linalg.eigvals(A_obs) assert all(e.real < 0 for e in eigs) + + +class TestGainSchedule: + def test_multiple_operating_points(self) -> None: + """Gain schedule over 3 points with varying spring constant.""" + + def linearize_spring(point: dict) -> tuple: + k = point["k"] + # dx/dt = v, dv/dt = -k*x - v + u (spring-mass-damper) + A = [[0.0, 1.0], [-k, -1.0]] + B = [[0.0], [1.0]] + C = [[1.0, 0.0]] + D = [[0.0]] + return (A, B, C, D) + + points = [{"k": 1.0}, {"k": 4.0}, {"k": 9.0}] + Q = [[1.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + + results = gain_schedule(linearize_spring, points, Q, R) + + assert len(results) == 3 + for point, K, E in results: + assert "k" in point + assert len(K) == 1 + assert len(K[0]) == 2 + assert all(e.real < 0 for e in E) + + # Gains should differ across operating points + K_values = [r[1][0][0] for r in results] + assert K_values[0] != pytest.approx(K_values[2], abs=0.01) diff --git a/packages/gds-domains/tests/test_transfer.py b/packages/gds-domains/tests/test_transfer.py index 63dd8d2..f4278c7 100644 --- a/packages/gds-domains/tests/test_transfer.py +++ b/packages/gds-domains/tests/test_transfer.py @@ -17,7 +17,6 @@ zeros, ) - # --------------------------------------------------------------------------- # Fixtures: canonical test systems # --------------------------------------------------------------------------- @@ -181,9 +180,7 @@ def test_first_order_pole(self, first_order_system: LinearizedSystem) -> None: assert p[0].real == pytest.approx(-1.0, abs=1e-6) assert abs(p[0].imag) < 1e-6 - def test_double_integrator_poles( - self, double_integrator: LinearizedSystem - ) -> None: + def test_double_integrator_poles(self, double_integrator: LinearizedSystem) -> None: tf_mat = ss_to_tf(double_integrator) tf = tf_mat.elements[0][0] p = poles(tf) @@ -280,7 +277,7 @@ def test_s_plus_t_equals_one(self) -> None: # S and T should have the same denominator assert len(s_tf.den) == len(t_tf.den) - for a, b in zip(s_tf.den, t_tf.den): + for a, b in zip(s_tf.den, t_tf.den, strict=False): assert a == pytest.approx(b, abs=1e-10) # S_num + T_num should equal the denominator @@ -290,8 +287,8 @@ def test_s_plus_t_equals_one(self) -> None: t_padded = [0.0] * (max_len - len(t_tf.num)) + t_tf.num d_padded = [0.0] * (max_len - len(s_tf.den)) + s_tf.den - summed = [a + b for a, b in zip(s_padded, t_padded)] - for a, b in zip(summed, d_padded): + summed = [a + b for a, b in zip(s_padded, t_padded, strict=False)] + for a, b in zip(summed, d_padded, strict=False): assert a == pytest.approx(b, abs=1e-10) def test_returns_all_six(self) -> None: @@ -300,3 +297,111 @@ def test_returns_all_six(self) -> None: gang = sensitivity(plant, controller) assert set(gang.keys()) == {"S", "T", "CS", "PS", "KS", "KPS"} + + def test_sensitivity_dc_values(self) -> None: + """Verify S(0) and T(0) at DC for P=1/(s+1), K=2. + + L(0) = P(0)*K(0) = 1*2 = 2 + S(0) = 1/(1+L(0)) = 1/3 + T(0) = L(0)/(1+L(0)) = 2/3 + """ + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) + controller = TransferFunction(num=[2.0], den=[1.0]) + + gang = sensitivity(plant, controller) + + # Evaluate S at s=0: ratio of constant terms + s_tf = gang["S"] + s_dc = s_tf.num[-1] / s_tf.den[-1] + assert s_dc == pytest.approx(1.0 / 3.0, abs=1e-10) + + t_tf = gang["T"] + t_dc = t_tf.num[-1] / t_tf.den[-1] + assert t_dc == pytest.approx(2.0 / 3.0, abs=1e-10) + + def test_cs_is_controller_times_s(self) -> None: + """CS = K*S: verify CS(0) = K(0)*S(0) = 2*(1/3) = 2/3.""" + plant = TransferFunction(num=[1.0], den=[1.0, 1.0]) + controller = TransferFunction(num=[2.0], den=[1.0]) + + gang = sensitivity(plant, controller) + cs_tf = gang["CS"] + cs_dc = cs_tf.num[-1] / cs_tf.den[-1] + assert cs_dc == pytest.approx(2.0 / 3.0, abs=1e-10) + + +# --------------------------------------------------------------------------- +# MIMO transfer functions +# --------------------------------------------------------------------------- + + +class TestMIMO: + def test_2x2_system(self) -> None: + """2-input, 2-output system: verify matrix dimensions.""" + ls = LinearizedSystem( + A=[[-1.0, 0.0], [0.0, -2.0]], + B=[[1.0, 0.0], [0.0, 1.0]], + C=[[1.0, 0.0], [0.0, 1.0]], + D=[[0.0, 0.0], [0.0, 0.0]], + x0=[0.0, 0.0], + u0=[0.0, 0.0], + state_names=["x1", "x2"], + input_names=["u1", "u2"], + output_names=["y1", "y2"], + ) + tf_mat = ss_to_tf(ls) + + # Should be 2x2 + assert len(tf_mat.elements) == 2 + assert len(tf_mat.elements[0]) == 2 + assert len(tf_mat.elements[1]) == 2 + + def test_2x2_diagonal_system_dc_gains(self) -> None: + """Diagonal A, identity B and C: verify DC gains. + + H11(s) = (s+2)/((s+1)(s+2)) → H11(0) = 2/2 = 1.0 + H22(s) = (s+1)/((s+1)(s+2)) → H22(0) = 1/2 = 0.5 + Note: denominator is the full characteristic polynomial for all elements. + """ + ls = LinearizedSystem( + A=[[-1.0, 0.0], [0.0, -2.0]], + B=[[1.0, 0.0], [0.0, 1.0]], + C=[[1.0, 0.0], [0.0, 1.0]], + D=[[0.0, 0.0], [0.0, 0.0]], + x0=[0.0, 0.0], + u0=[0.0, 0.0], + state_names=["x1", "x2"], + input_names=["u1", "u2"], + output_names=["y1", "y2"], + ) + tf_mat = ss_to_tf(ls) + + # DC gain = num(0)/den(0) = constant_term_num / constant_term_den + h11 = tf_mat.elements[0][0] + h11_dc = h11.num[-1] / h11.den[-1] + assert h11_dc == pytest.approx(1.0, abs=1e-8) + + h22 = tf_mat.elements[1][1] + h22_dc = h22.num[-1] / h22.den[-1] + assert h22_dc == pytest.approx(0.5, abs=1e-8) + + def test_nonzero_D_feedthrough(self) -> None: + """System with D != 0: verify DC gain includes feedthrough.""" + ls = LinearizedSystem( + A=[[-1.0]], + B=[[1.0]], + C=[[1.0]], + D=[[2.0]], + x0=[0.0], + u0=[0.0], + state_names=["x"], + input_names=["u"], + output_names=["y"], + ) + tf_mat = ss_to_tf(ls) + tf = tf_mat.elements[0][0] + + # H(s) = 1/(s+1) + 2 = (2s + 3)/(s + 1) + # DC gain H(0) = 3/1 = 3 + dc_gain = tf.num[-1] / tf.den[-1] + assert dc_gain == pytest.approx(3.0, abs=1e-8) diff --git a/packages/gds-proof/CLAUDE.md b/packages/gds-proof/CLAUDE.md index 49befec..918e824 100644 --- a/packages/gds-proof/CLAUDE.md +++ b/packages/gds-proof/CLAUDE.md @@ -161,6 +161,24 @@ model.canonical_dict() --- +## Lyapunov Stability Proofs (`analysis/lyapunov.py`) + +Control-theoretic stability templates using SymPy. Does not use the five-strategy +prover from `symbolic.py` — performs its own Hessian-based definiteness checks. + +| Function | What it proves | +|----------|---------------| +| `lyapunov_candidate(V, f, states)` | V(x) > 0 and dV/dt < 0 (continuous) or ΔV < 0 (discrete) | +| `quadratic_lyapunov(P, A)` | V = x'Px via P > 0 and A'P + PA < 0 (eigenvalue check) | +| `find_quadratic_lyapunov(A, Q)` | Solve A'P + PA = -Q for P (SymPy linear system) | +| `passivity_certificate(V, s, f, h)` | dV/dt ≤ s(u, y) for dissipativity | + +Results are `LyapunovResult` / `PassivityResult` with `Literal["PROVED", "FAILED", "INCONCLUSIVE"]` status fields. + +**Limitation:** Only handles quadratic forms (constant Hessian) definitively. Non-quadratic V returns INCONCLUSIVE. For systems with n > 4, `find_quadratic_lyapunov` may be slow (use `scipy.linalg.solve_continuous_lyapunov` in `gds-analysis` instead). + +--- + ## Commands ```bash diff --git a/packages/gds-proof/gds_proof/analysis/lyapunov.py b/packages/gds-proof/gds_proof/analysis/lyapunov.py index 9cd52a8..04afe4c 100644 --- a/packages/gds-proof/gds_proof/analysis/lyapunov.py +++ b/packages/gds-proof/gds_proof/analysis/lyapunov.py @@ -9,8 +9,8 @@ For discrete-time systems: V(x) > 0 and V(f(x)) - V(x) < 0 -All proofs delegate to the existing five-strategy SymPy simplification -engine in ``gds_proof.analysis.symbolic``. +Uses SymPy's simplification engine directly for positive-definiteness +and decrease condition checks via Hessian eigenvalue analysis. """ from __future__ import annotations diff --git a/packages/gds-symbolic/CLAUDE.md b/packages/gds-symbolic/CLAUDE.md index de2302c..3341a8f 100644 --- a/packages/gds-symbolic/CLAUDE.md +++ b/packages/gds-symbolic/CLAUDE.md @@ -16,6 +16,8 @@ | `model.py` | `SymbolicControlModel(ControlModel)` — adds symbolic equations + validation | | `compile.py` | `compile_to_ode(model)` → `(ODEFunction, state_order)` via `parse_expr` + `lambdify` | | `linearize.py` | `linearize(model, x0, u0)` → `LinearizedSystem(A, B, C, D)` via Jacobians | +| `transfer.py` | `ss_to_tf()`, `poles()`, `zeros()`, `controllability_matrix()`, `sensitivity()` (Gang of Six) | +| `delay.py` | `pade_approximation()`, `delay_system()` — Padé time delay modeling | | `errors.py` | `SymbolicError(CSError)` — inherits control domain error hierarchy | ### Security @@ -28,10 +30,25 @@ Expression parsing uses `sympy.parsing.sympy_parser.parse_expr` with a restricte - `sympy.Expr` (parsed) → R2, transient - Lambdified callable → R3, opaque +### Transfer function analysis (`transfer.py`) + +Builds on `LinearizedSystem` (A, B, C, D as `list[list[float]]`): +- `ss_to_tf(ls)` → `TransferFunctionMatrix` via adjugate/determinant method +- `poles(tf)`, `zeros(tf)` → `list[complex]` via SymPy roots +- `controllability_matrix(ls)`, `observability_matrix(ls)` → Kalman rank matrices +- `sensitivity(plant, controller)` → Gang of Six: S, T, CS, PS, KS, KPS + +### Padé time delay (`delay.py`) + +- `pade_approximation(delay, order)` → all-pass `TransferFunction` approximating e^{-sτ} +- `delay_system(tf, delay, order)` → cascaded TF with Padé delay + ### Known limitations - Inputs resolved from constant `params` — no time-varying signals (see gds-continuous CLAUDE.md) - `compile()` and `compile_system()` remain structural — symbolic equations don't appear in GDSSpec +- `poles()`/`zeros()` use SymPy's symbolic root finder — may return incomplete results for degree > 4 +- Transfer functions use the full characteristic polynomial as denominator (MIMO elements share the same denominator) ## Commands diff --git a/packages/gds-viz/CLAUDE.md b/packages/gds-viz/CLAUDE.md index 5b29ccc..7523acf 100644 --- a/packages/gds-viz/CLAUDE.md +++ b/packages/gds-viz/CLAUDE.md @@ -6,7 +6,7 @@ - **Import**: `import gds_viz` - **Dependencies**: `gds-framework>=0.2.3` -- **Optional**: `[phase]` for matplotlib + numpy + gds-continuous (phase portraits) +- **Optional**: `[phase]` for matplotlib + numpy + gds-continuous (phase portraits), `[control]` for matplotlib + numpy (frequency/response plots), `[plots]` for all plotting ## Architecture @@ -33,6 +33,25 @@ Supports >2D systems via `fixed_states` projection (e.g., Lorenz z=25 slice). +### Frequency response plots (`gds_viz.frequency`, requires `[control]`) + +| Function | Purpose | +|----------|---------| +| `bode_plot(omega, mag_db, phase_deg)` | Magnitude + phase with margin annotations | +| `nyquist_plot(real, imag)` | Complex plane with unit circle and critical point | +| `nichols_plot(phase_deg, mag_db)` | Open-loop phase vs. gain with M-circles | +| `root_locus_plot(num, den)` | Pole migration with gain sweep | + +### Time-domain response plots (`gds_viz.response`, requires `[control]`) + +| Function | Purpose | +|----------|---------| +| `step_response_plot(times, values, metrics=...)` | Annotated step response with metric overlays | +| `impulse_response_plot(times, values)` | Basic impulse response | +| `compare_responses(responses)` | Multi-response overlay comparison | + +All plot functions accept plain `list[float]` data — no coupling to gds-analysis types. Lazy-imported via `__getattr__`. + ## Commands ```bash diff --git a/packages/gds-viz/gds_viz/frequency.py b/packages/gds-viz/gds_viz/frequency.py index 3dc5daf..0250499 100644 --- a/packages/gds-viz/gds_viz/frequency.py +++ b/packages/gds-viz/gds_viz/frequency.py @@ -208,7 +208,6 @@ def nichols_plot( *, title: str = "Nichols Chart", m_circles: bool = True, - n_circles: bool = False, ax: Any | None = None, figsize: tuple[float, float] = (10, 8), ) -> Any: @@ -224,8 +223,6 @@ def nichols_plot( Plot title. m_circles : bool If True, draw M-circles (constant closed-loop magnitude contours). - n_circles : bool - If True, draw N-circles (constant closed-loop phase contours). ax : Axes | None Existing axes. figsize : tuple[float, float] diff --git a/packages/gds-viz/gds_viz/response.py b/packages/gds-viz/gds_viz/response.py index fcac33c..ddfb750 100644 --- a/packages/gds-viz/gds_viz/response.py +++ b/packages/gds-viz/gds_viz/response.py @@ -167,7 +167,11 @@ def _annotate_metrics(ax: Any, t: Any, y: Any, metrics: Any) -> None: def _find_crossing_time(t: Any, y: Any, threshold: float) -> float: - """Find first time where y crosses threshold.""" + """Find first time where y crosses threshold via linear interpolation. + + Note: identical logic exists in ``gds_analysis.response._interpolate_crossing``. + Duplicated here to avoid a hard dependency on gds-analysis from gds-viz. + """ for i in range(1, len(y)): if (y[i - 1] <= threshold <= y[i]) or (y[i - 1] >= threshold >= y[i]): dv = y[i] - y[i - 1] diff --git a/packages/gds-viz/tests/test_frequency.py b/packages/gds-viz/tests/test_frequency.py new file mode 100644 index 0000000..5db3826 --- /dev/null +++ b/packages/gds-viz/tests/test_frequency.py @@ -0,0 +1,129 @@ +"""Tests for frequency response and root locus plots.""" + +import pytest + +from gds_viz.frequency import ( + bode_plot, + nichols_plot, + nyquist_plot, + root_locus_plot, +) + + +@pytest.fixture(autouse=True) +def _mpl_backend(): + """Use non-interactive backend for CI.""" + import matplotlib + + matplotlib.use("Agg") + + +class TestBodePlot: + def test_returns_figure_with_two_axes(self) -> None: + omega = [0.1, 1.0, 10.0, 100.0] + mag_db = [0.0, -3.0, -20.0, -40.0] + phase_deg = [0.0, -45.0, -90.0, -90.0] + + fig = bode_plot(omega, mag_db, phase_deg) + + assert fig is not None + axes = fig.get_axes() + assert len(axes) == 2 + + def test_magnitude_axis_is_semilog(self) -> None: + fig = bode_plot([0.1, 10.0], [-3.0, -20.0], [-45.0, -90.0]) + ax_mag = fig.get_axes()[0] + assert ax_mag.get_xscale() == "log" + + def test_with_margin_annotations(self) -> None: + fig = bode_plot( + [0.1, 1.0, 10.0], + [20.0, 0.0, -20.0], + [-90.0, -135.0, -180.0], + gain_margin=(6.0, 5.0), + phase_margin=(45.0, 1.0), + ) + assert fig is not None + + def test_with_existing_axes(self) -> None: + import matplotlib.pyplot as plt + + fig, (ax1, ax2) = plt.subplots(2, 1) + result = bode_plot([0.1, 10.0], [-3.0, -20.0], [-45.0, -90.0], ax=(ax1, ax2)) + assert result is fig + + def test_custom_title(self) -> None: + fig = bode_plot([1.0], [0.0], [-90.0], title="My Custom Bode") + ax = fig.get_axes()[0] + assert ax.get_title() == "My Custom Bode" + + +class TestNyquistPlot: + def test_returns_figure(self) -> None: + real = [1.0, 0.5, 0.0, -0.5] + imag = [0.0, -0.5, -1.0, -0.5] + + fig = nyquist_plot(real, imag) + assert fig is not None + + def test_has_critical_point_marker(self) -> None: + fig = nyquist_plot([1.0, 0.0], [0.0, -1.0]) + ax = fig.get_axes()[0] + # Should have line artists for the curve + mirror + unit circle + critical point + assert len(ax.lines) >= 3 + + def test_equal_aspect(self) -> None: + fig = nyquist_plot([1.0, 0.0], [0.0, -1.0]) + ax = fig.get_axes()[0] + # matplotlib may return "equal", "equalxy", or 1.0 depending on version + assert ax.get_aspect() in ("equal", "equalxy", 1.0) + + +class TestNicholsPlot: + def test_returns_figure(self) -> None: + phase = [-90.0, -135.0, -180.0, -225.0] + mag = [20.0, 6.0, 0.0, -6.0] + + fig = nichols_plot(phase, mag) + assert fig is not None + + def test_x_axis_limits(self) -> None: + fig = nichols_plot([-100.0, -200.0], [10.0, -10.0]) + ax = fig.get_axes()[0] + xlim = ax.get_xlim() + assert xlim[0] == -360 + assert xlim[1] == 0 + + def test_without_m_circles(self) -> None: + fig = nichols_plot([-90.0, -180.0], [0.0, -10.0], m_circles=False) + assert fig is not None + + +class TestRootLocusPlot: + def test_returns_figure(self) -> None: + # 1/(s^2 + s + 1): two poles + num = [1.0] + den = [1.0, 1.0, 1.0] + + fig = root_locus_plot(num, den) + assert fig is not None + + def test_marks_poles_and_zeros(self) -> None: + # (s+2)/(s^2+3s+2) = (s+2)/((s+1)(s+2)): 2 poles, 1 zero + num = [1.0, 2.0] + den = [1.0, 3.0, 2.0] + + fig = root_locus_plot(num, den, mark_poles=True, mark_zeros=True) + ax = fig.get_axes()[0] + # Should have line artists for trajectories + poles marker + zeros marker + assert len(ax.lines) >= 3 + + def test_custom_gains(self) -> None: + fig = root_locus_plot([1.0], [1.0, 2.0, 1.0], gains=[0.0, 0.5, 1.0, 5.0, 10.0]) + assert fig is not None + + def test_equal_aspect(self) -> None: + fig = root_locus_plot([1.0], [1.0, 1.0]) + ax = fig.get_axes()[0] + # matplotlib may return "equal", "equalxy", or 1.0 depending on version + assert ax.get_aspect() in ("equal", "equalxy", 1.0) diff --git a/packages/gds-viz/tests/test_response_viz.py b/packages/gds-viz/tests/test_response_viz.py new file mode 100644 index 0000000..020acaa --- /dev/null +++ b/packages/gds-viz/tests/test_response_viz.py @@ -0,0 +1,140 @@ +"""Tests for step/impulse response visualization.""" + +import math + +import pytest + +from gds_viz.response import ( + compare_responses, + impulse_response_plot, + step_response_plot, +) + + +@pytest.fixture(autouse=True) +def _mpl_backend(): + """Use non-interactive backend for CI.""" + import matplotlib + + matplotlib.use("Agg") + + +@pytest.fixture() +def sample_step_data() -> tuple[list[float], list[float]]: + """First-order step response: y = 1 - e^{-t}.""" + n = 200 + t_max = 8.0 + times = [i * t_max / (n - 1) for i in range(n)] + values = [1.0 - math.exp(-t) for t in times] + return times, values + + +class TestStepResponsePlot: + def test_returns_figure(self, sample_step_data: tuple) -> None: + times, values = sample_step_data + fig = step_response_plot(times, values) + assert fig is not None + + def test_with_setpoint(self, sample_step_data: tuple) -> None: + times, values = sample_step_data + fig = step_response_plot(times, values, setpoint=1.0) + ax = fig.get_axes()[0] + # Should have response line + setpoint line + assert len(ax.lines) >= 2 + + def test_with_metrics_annotation(self, sample_step_data: tuple) -> None: + """Verify metrics annotations don't crash and add artists.""" + from gds_analysis.response import StepMetrics + + times, values = sample_step_data + metrics = StepMetrics( + rise_time=2.2, + settling_time=4.0, + overshoot_pct=0.0, + peak_time=8.0, + peak_value=0.999, + steady_state_value=1.0, + steady_state_error=0.0, + ) + fig = step_response_plot(times, values, setpoint=1.0, metrics=metrics) + assert fig is not None + # Should have patches (settling band, rise time span) + ax = fig.get_axes()[0] + assert len(ax.patches) >= 1 # At least the settling band + + def test_with_overshoot_metrics(self) -> None: + """Verify overshoot annotation renders.""" + from gds_analysis.response import StepMetrics + + times = [0.0, 0.5, 1.0, 1.5, 2.0, 5.0] + values = [0.0, 0.8, 1.2, 1.05, 1.0, 1.0] + metrics = StepMetrics( + rise_time=0.6, + settling_time=2.0, + overshoot_pct=20.0, + peak_time=1.0, + peak_value=1.2, + steady_state_value=1.0, + steady_state_error=0.0, + ) + fig = step_response_plot(times, values, metrics=metrics) + assert fig is not None + + def test_custom_title(self, sample_step_data: tuple) -> None: + times, values = sample_step_data + fig = step_response_plot(times, values, title="My Step") + ax = fig.get_axes()[0] + assert ax.get_title() == "My Step" + + def test_existing_axes(self, sample_step_data: tuple) -> None: + import matplotlib.pyplot as plt + + times, values = sample_step_data + fig, ax = plt.subplots() + result = step_response_plot(times, values, ax=ax) + assert result is fig + + +class TestImpulseResponsePlot: + def test_returns_figure(self) -> None: + times = [0.0, 1.0, 2.0, 3.0, 4.0] + values = [1.0, 0.37, 0.14, 0.05, 0.02] + fig = impulse_response_plot(times, values) + assert fig is not None + + def test_has_zero_line(self) -> None: + fig = impulse_response_plot([0.0, 1.0], [1.0, 0.0]) + ax = fig.get_axes()[0] + # Should have the response line + zero reference line + assert len(ax.lines) >= 2 + + +class TestCompareResponses: + def test_returns_figure_with_multiple_lines(self) -> None: + responses = [ + ([0.0, 1.0, 2.0], [0.0, 0.6, 0.9], "Slow"), + ([0.0, 1.0, 2.0], [0.0, 0.9, 1.0], "Fast"), + ([0.0, 1.0, 2.0], [0.0, 1.2, 1.0], "Overshoot"), + ] + fig = compare_responses(responses) + ax = fig.get_axes()[0] + assert len(ax.lines) == 3 + + def test_single_response(self) -> None: + fig = compare_responses( + [ + ([0.0, 1.0], [0.0, 1.0], "Only one"), + ] + ) + assert fig is not None + + def test_legend_present(self) -> None: + fig = compare_responses( + [ + ([0.0, 1.0], [0.0, 1.0], "A"), + ([0.0, 1.0], [0.0, 0.5], "B"), + ] + ) + ax = fig.get_axes()[0] + legend = ax.get_legend() + assert legend is not None From 74f5a7af3a8a8b1e7a95b60762a9162f06cf1b40 Mon Sep 17 00:00:00 2001 From: rohan Date: Tue, 7 Apr 2026 19:47:40 +0530 Subject: [PATCH 3/3] =?UTF-8?q?feat:=20add=20DC=20motor=20control=20tutori?= =?UTF-8?q?al=20=E2=80=94=208=20lessons,=20beginner=20to=20advanced?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Progressive tutorial using a DC motor position controller to teach classical control theory concepts with GDS. Each lesson is standalone (repeats motor parameters), builds on the same physical system, and exercises the new control theory analysis stack from #198. Lessons: 1. lesson1_plant — SymbolicControlModel, ODE simulation, phase portrait 2. lesson2_p_control — P controller, StepMetrics, gain sweep 3. lesson3_pid_transfer — PID, linearize, transfer functions, Bode 4. lesson4_disturbance — feedforward, Gang of Six sensitivity analysis 5. lesson5_delay — Pade approximation, gain/phase margins, Nyquist 6. lesson6_lqr — LQR vs PID, Kalman filter, gain scheduling 7. lesson7_discrete — Tustin/ZOH discretization, gds-sim execution 8. lesson8_lyapunov — Lyapunov proof, passivity certificate 31 tests in test_control_tutorial.py, 884 total gds-examples tests pass. --- .../gds_examples/control_tutorial/__init__.py | 16 + .../control_tutorial/lesson1_plant.py | 158 +++++++++ .../control_tutorial/lesson2_p_control.py | 146 ++++++++ .../control_tutorial/lesson3_pid_transfer.py | 229 ++++++++++++ .../control_tutorial/lesson4_disturbance.py | 160 +++++++++ .../control_tutorial/lesson5_delay.py | 164 +++++++++ .../control_tutorial/lesson6_lqr.py | 199 +++++++++++ .../control_tutorial/lesson7_discrete.py | 220 ++++++++++++ .../control_tutorial/lesson8_lyapunov.py | 182 ++++++++++ .../tests/test_control_tutorial.py | 331 ++++++++++++++++++ 10 files changed, 1805 insertions(+) create mode 100644 packages/gds-examples/gds_examples/control_tutorial/__init__.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson1_plant.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson2_p_control.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson3_pid_transfer.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson4_disturbance.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson5_delay.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson6_lqr.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson7_discrete.py create mode 100644 packages/gds-examples/gds_examples/control_tutorial/lesson8_lyapunov.py create mode 100644 packages/gds-examples/tests/test_control_tutorial.py diff --git a/packages/gds-examples/gds_examples/control_tutorial/__init__.py b/packages/gds-examples/gds_examples/control_tutorial/__init__.py new file mode 100644 index 0000000..4dbcec7 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/__init__.py @@ -0,0 +1,16 @@ +"""Classical Control Theory with GDS -- progressive tutorial. + +A DC motor position controller serves as the running example across +eight lessons, each introducing new control theory concepts and the +GDS modules that implement them. + +Lessons: + 1. lesson1_plant -- Model the DC motor (SymbolicControlModel, ODE sim) + 2. lesson2_p_control -- Proportional control and step response metrics + 3. lesson3_pid_transfer -- PID, transfer functions, Bode analysis + 4. lesson4_disturbance -- Disturbance rejection and sensitivity (Gang of Six) + 5. lesson5_delay -- Sensor delay, Pade approximation, stability margins + 6. lesson6_lqr -- LQR optimal control, Kalman filter, gain scheduling + 7. lesson7_discrete -- Discretization and discrete-time simulation (gds-sim) + 8. lesson8_lyapunov -- Lyapunov stability proof and passivity certificate +""" diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson1_plant.py b/packages/gds-examples/gds_examples/control_tutorial/lesson1_plant.py new file mode 100644 index 0000000..adc936a --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson1_plant.py @@ -0,0 +1,158 @@ +"""Lesson 1 -- Modeling a DC Motor Plant. + +Build a symbolic model of a DC motor that controls angular position. +The motor converts voltage to torque, which spins a shaft against +viscous friction. We declare the structural GDS spec (blocks, wiring, +roles) and the behavioral ODEs in a single SymbolicControlModel, then +simulate the open-loop response to a constant voltage step. + +New concepts: + - SymbolicControlModel (from gds_domains.symbolic) + - StateEquation / OutputEquation for symbolic ODEs + - compile_model() / compile_to_system() for GDS compilation + - to_ode_function() -> ODEFunction for simulation + - ODEModel / ODESimulation (from gds_continuous) + +GDS Decomposition: + X = (theta, omega) -- angular position, angular velocity + U = V -- voltage input + g = voltage_command -- trivial pass-through (open-loop) + f = motor_dynamics -- dtheta/dt = omega, domega/dt = ... + Theta = {J, b, K_t} -- inertia, damping, torque constant + +Composition: + voltage_input >> voltage_command >> motor_dynamics +""" + +from gds.ir.models import SystemIR +from gds.spec import GDSSpec +from gds_domains.control.dsl.compile import compile_model, compile_to_system +from gds_domains.control.dsl.elements import Controller, Input, Sensor, State +from gds_domains.symbolic.elements import OutputEquation, StateEquation +from gds_domains.symbolic.model import SymbolicControlModel + +# --------------------------------------------------------------------------- +# Motor parameters (repeated in every lesson for standalone use) +# --------------------------------------------------------------------------- + +J = 0.01 # moment of inertia (kg*m^2) +b = 0.1 # viscous damping coefficient (N*m*s/rad) +K_t = 0.01 # motor torque constant (N*m/A, simplified V->torque) + + +# --------------------------------------------------------------------------- +# Build functions +# --------------------------------------------------------------------------- + + +def build_model() -> SymbolicControlModel: + """Declare the DC motor as a SymbolicControlModel. + + Structural layer (what GDS sees): + 2 states, 1 input, 1 sensor, 1 controller + + Behavioral layer (what the ODE integrator sees): + dtheta/dt = omega + domega/dt = -(b/J)*omega + (K_t/J)*V + """ + return SymbolicControlModel( + name="DC Motor Plant", + states=[ + State(name="theta", initial=0.0), + State(name="omega", initial=0.0), + ], + inputs=[ + Input(name="V"), + ], + sensors=[ + Sensor(name="position_sensor", observes=["theta"]), + ], + controllers=[ + Controller( + name="voltage_command", + reads=["position_sensor", "V"], + drives=["theta", "omega"], + ), + ], + state_equations=[ + StateEquation(state_name="theta", expr_str="omega"), + StateEquation( + state_name="omega", + expr_str=f"-({b}/{J})*omega + ({K_t}/{J})*V", + ), + ], + output_equations=[ + OutputEquation(sensor_name="position_sensor", expr_str="theta"), + ], + symbolic_params=[], + description="DC motor angular position control plant", + ) + + +def build_spec() -> GDSSpec: + """Compile the SymbolicControlModel to a full GDSSpec.""" + return compile_model(build_model()) + + +def build_system() -> SystemIR: + """Compile to SystemIR with automatic composition tree.""" + return compile_to_system(build_model()) + + +def simulate_open_loop( + voltage: float = 1.0, + t_end: float = 2.0, +) -> tuple[list[float], list[float], list[float]]: + """Simulate open-loop step response (constant voltage). + + Returns (times, theta_values, omega_values). + """ + from gds_continuous import ODEModel, ODESimulation + + model = build_model() + ode_fn, state_order = model.to_ode_function() + + ode_model = ODEModel( + state_names=state_order, + initial_state={"theta": 0.0, "omega": 0.0}, + rhs=ode_fn, + params={"V": [voltage]}, + ) + sim = ODESimulation( + model=ode_model, + t_span=(0.0, t_end), + solver="RK45", + ) + results = sim.run() + + return ( + results.times, + results.state_array("theta"), + results.state_array("omega"), + ) + + +if __name__ == "__main__": + # --- Structural layer --- + model = build_model() + print("=== Lesson 1: DC Motor Plant ===\n") + print(f"Model: {model.name}") + print(f" States: {[s.name for s in model.states]}") + print(f" Inputs: {[i.name for i in model.inputs]}") + print(f" Sensors: {[s.name for s in model.sensors]}") + print(f" Controllers: {[c.name for c in model.controllers]}") + + spec = build_spec() + print(f"\nGDSSpec: {spec.name}") + print(f" Blocks: {list(spec.blocks.keys())}") + print(f" Errors: {spec.validate_spec()}") + + system = build_system() + print(f"\nSystemIR: {len(system.blocks)} blocks, {len(system.wirings)} wirings") + + # --- Behavioral layer --- + print("\n--- Open-loop simulation (V=1.0, 2 seconds) ---") + times, theta, omega = simulate_open_loop(voltage=1.0, t_end=2.0) + print(f" Final theta = {theta[-1]:.4f} rad") + print(f" Final omega = {omega[-1]:.4f} rad/s") + print(" (Motor spins up to steady velocity, position ramps)") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson2_p_control.py b/packages/gds-examples/gds_examples/control_tutorial/lesson2_p_control.py new file mode 100644 index 0000000..33b9e2d --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson2_p_control.py @@ -0,0 +1,146 @@ +"""Lesson 2 -- Proportional Control and Performance Metrics. + +Close the loop with a proportional (P) controller: V = Kp * (theta_ref - theta). +Simulate the step response for different gain values and extract standard +performance metrics: rise time, settling time, overshoot, steady-state error. + +New concepts: + - Closed-loop control via ODE parameter (Kp, theta_ref) + - step_response_metrics() / StepMetrics (from gds_analysis.response) + - step_response_plot() with metric annotations (from gds_viz.response) + - compare_responses() for gain sweep visualization + +GDS Decomposition: + X = (theta, omega) + U = theta_ref (reference position) + g = p_controller: V = Kp * (theta_ref - theta) + f = motor_dynamics + Theta = {J, b, K_t, Kp, theta_ref} + +Composition: + (reference | position_sensor) >> p_controller >> motor_dynamics + .loop(dynamics -> position_sensor) +""" + +from typing import Any + +from gds_domains.control.dsl.elements import Controller, Input, Sensor, State +from gds_domains.symbolic.elements import OutputEquation, StateEquation +from gds_domains.symbolic.model import SymbolicControlModel + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Build functions +# --------------------------------------------------------------------------- + + +def build_model() -> SymbolicControlModel: + """DC motor with proportional position controller. + + The controller law V = Kp * (theta_ref - theta) is embedded in the + state equation for omega. In GDS terms, the controller is a Policy + block whose behavioral implementation is this proportional law. + """ + return SymbolicControlModel( + name="DC Motor P Control", + states=[ + State(name="theta", initial=0.0), + State(name="omega", initial=0.0), + ], + inputs=[ + Input(name="theta_ref"), + ], + sensors=[ + Sensor(name="position_sensor", observes=["theta"]), + ], + controllers=[ + Controller( + name="p_controller", + reads=["position_sensor", "theta_ref"], + drives=["theta", "omega"], + ), + ], + state_equations=[ + StateEquation(state_name="theta", expr_str="omega"), + StateEquation( + state_name="omega", + expr_str=(f"-({b}/{J})*omega + ({K_t}/{J})*Kp*(theta_ref - theta)"), + ), + ], + output_equations=[ + OutputEquation(sensor_name="position_sensor", expr_str="theta"), + ], + symbolic_params=["Kp"], + description="DC motor with proportional position control", + ) + + +def simulate_p_control( + kp: float = 5.0, + theta_ref: float = 1.0, + t_end: float = 5.0, +) -> tuple[list[float], list[float]]: + """Simulate closed-loop P control step response. + + Returns (times, theta_values). + """ + from gds_continuous import ODEModel, ODESimulation + + model = build_model() + ode_fn, state_order = model.to_ode_function() + + ode_model = ODEModel( + state_names=state_order, + initial_state={"theta": 0.0, "omega": 0.0}, + rhs=ode_fn, + params={"theta_ref": [theta_ref], "Kp": [kp]}, + ) + sim = ODESimulation(model=ode_model, t_span=(0.0, t_end), solver="RK45") + results = sim.run() + return results.times, results.state_array("theta") + + +def analyze_step_response( + kp: float = 5.0, + theta_ref: float = 1.0, +) -> tuple[list[float], list[float], Any]: + """Simulate and compute step response metrics. + + Returns (times, theta_values, StepMetrics). + """ + from gds_analysis.response import step_response_metrics + + times, theta = simulate_p_control(kp=kp, theta_ref=theta_ref) + metrics = step_response_metrics(times, theta, setpoint=theta_ref) + return times, theta, metrics + + +if __name__ == "__main__": + print("=== Lesson 2: Proportional Control ===\n") + + # --- Gain sweep --- + gains = [1.0, 5.0, 10.0] + print(f"{'Kp':>6} {'Rise':>8} {'Settle':>8} {'OS%':>8} {'SSE':>8}") + print("-" * 50) + + responses = [] + for kp in gains: + times, theta, metrics = analyze_step_response(kp=kp) + print( + f"{kp:6.1f} {metrics.rise_time:8.3f} " + f"{metrics.settling_time:8.3f} {metrics.overshoot_pct:8.2f} " + f"{metrics.steady_state_error:8.4f}" + ) + responses.append((times, theta, f"Kp={kp}")) + + print("\n Higher Kp -> faster rise but more overshoot.") + print(" P control alone has steady-state error for this 2nd-order plant.") + print(" (Lesson 3 adds integral action to eliminate SSE.)") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson3_pid_transfer.py b/packages/gds-examples/gds_examples/control_tutorial/lesson3_pid_transfer.py new file mode 100644 index 0000000..6c6fd74 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson3_pid_transfer.py @@ -0,0 +1,229 @@ +"""Lesson 3 -- PID, Transfer Functions, and Bode Analysis. + +Upgrade from P to PID control, then cross into the frequency domain. +Linearize the plant to get (A, B, C, D) matrices, convert to a transfer +function, and analyze poles, zeros, controllability, and the Bode plot. + +New concepts: + - linearize() -> LinearizedSystem (A, B, C, D) + - ss_to_tf() -> TransferFunctionMatrix + - poles(), zeros(), is_controllable(), is_observable() + - frequency_response() for numerical Bode data + - bode_plot() for visualization + +GDS Decomposition: + X = (theta, omega) + U = theta_ref + g = pid_controller: V = Kp*e + Ki*integral(e) + Kd*de/dt + f = motor_dynamics + Theta = {J, b, K_t, Kp, Ki, Kd} + +Composition: + (reference | position_sensor) >> pid_controller >> motor_dynamics + .loop(dynamics -> position_sensor) +""" + +from gds_domains.control.dsl.elements import Controller, Input, Sensor, State +from gds_domains.symbolic.elements import OutputEquation, StateEquation +from gds_domains.symbolic.model import SymbolicControlModel + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Build functions +# --------------------------------------------------------------------------- + + +def build_model() -> SymbolicControlModel: + """DC motor with PID controller embedded in the ODE. + + The integral of the error is modeled as an augmented state + variable 'error_int'. This lets the PID law appear in the + state equations without requiring a separate transfer function. + """ + # PID: V = Kp*(ref-theta) + Ki*error_int + Kd*(0 - omega) + # d(error_int)/dt = theta_ref - theta + # Note: derivative term uses -omega (derivative of error = -omega + # when theta_ref is constant) + return SymbolicControlModel( + name="DC Motor PID", + states=[ + State(name="theta", initial=0.0), + State(name="omega", initial=0.0), + ], + inputs=[ + Input(name="theta_ref"), + ], + sensors=[ + Sensor(name="position_sensor", observes=["theta"]), + ], + controllers=[ + Controller( + name="pid_controller", + reads=["position_sensor", "theta_ref"], + drives=["theta", "omega"], + ), + ], + state_equations=[ + StateEquation(state_name="theta", expr_str="omega"), + StateEquation( + state_name="omega", + expr_str=( + f"-({b}/{J})*omega" + f" + ({K_t}/{J})*(Kp*(theta_ref - theta)" + f" + Ki*error_int" + f" + Kd*(0 - omega))" + ), + ), + ], + output_equations=[ + OutputEquation(sensor_name="position_sensor", expr_str="theta"), + ], + symbolic_params=["Kp", "Ki", "Kd", "error_int"], + description="DC motor with PID position control", + ) + + +def get_plant_linearization(): + """Linearize the open-loop plant (no controller) at the origin. + + Returns LinearizedSystem with A, B, C, D for the 2-state motor. + This is the PLANT only — no controller in the loop. + """ + from gds_domains.symbolic.linearize import LinearizedSystem + + # Analytical state-space for the DC motor plant: + # A = [[0, 1], [0, -b/J]] + # B = [[0], [K_t/J]] + # C = [[1, 0]] (position output) + # D = [[0]] + return LinearizedSystem( + A=[[0.0, 1.0], [0.0, -b / J]], + B=[[0.0], [K_t / J]], + C=[[1.0, 0.0]], + D=[[0.0]], + x0=[0.0, 0.0], + u0=[0.0], + state_names=["theta", "omega"], + input_names=["V"], + output_names=["theta_out"], + ) + + +def analyze_transfer_function(): + """Compute open-loop plant transfer function and properties. + + Returns (tf, plant_poles, plant_zeros, controllable, observable). + """ + from gds_domains.symbolic.transfer import ( + is_controllable, + is_observable, + poles, + ss_to_tf, + zeros, + ) + + ls = get_plant_linearization() + tf_mat = ss_to_tf(ls) + tf = tf_mat.elements[0][0] + + return ( + tf, + poles(tf), + zeros(tf), + is_controllable(ls), + is_observable(ls), + ) + + +def simulate_pid( + kp: float = 10.0, + ki: float = 5.0, + kd: float = 0.5, + theta_ref: float = 1.0, + t_end: float = 5.0, +) -> tuple[list[float], list[float]]: + """Simulate PID-controlled motor. + + Uses a manually constructed ODE that includes the error integral + as a third state variable (augmented system). + """ + from gds_continuous import ODEModel, ODESimulation + + def pid_ode(t, y, params): + theta, omega, err_int = y + ref = params.get("theta_ref", 1.0) + kp_ = params.get("Kp", 10.0) + ki_ = params.get("Ki", 5.0) + kd_ = params.get("Kd", 0.5) + + error = ref - theta + v = kp_ * error + ki_ * err_int + kd_ * (0 - omega) + + dtheta = omega + domega = -(b / J) * omega + (K_t / J) * v + derr_int = error + return [dtheta, domega, derr_int] + + ode_model = ODEModel( + state_names=["theta", "omega", "error_int"], + initial_state={"theta": 0.0, "omega": 0.0, "error_int": 0.0}, + rhs=pid_ode, + params={ + "theta_ref": [theta_ref], + "Kp": [kp], + "Ki": [ki], + "Kd": [kd], + }, + ) + sim = ODESimulation(model=ode_model, t_span=(0.0, t_end), solver="RK45") + results = sim.run() + return results.times, results.state_array("theta") + + +if __name__ == "__main__": + from gds_analysis.linear import eigenvalues, frequency_response, is_stable + from gds_analysis.response import step_response_metrics + + print("=== Lesson 3: PID, Transfer Functions, Bode ===\n") + + # --- Plant linearization --- + ls = get_plant_linearization() + print("Plant state-space (open-loop):") + print(f" A = {ls.A}") + print(f" B = {ls.B}") + print(f" C = {ls.C}") + print(f" D = {ls.D}") + + eigs = eigenvalues(ls.A) + print(f" Eigenvalues: {[f'{e.real:.2f}' for e in eigs]}") + print(f" Open-loop stable: {is_stable(ls.A)}") + + # --- Transfer function --- + tf, p, z, ctrl, obs = analyze_transfer_function() + print(f"\nPlant TF: num={tf.num}, den={tf.den}") + print(f" Poles: {[f'{pi.real:.2f}+{pi.imag:.2f}j' for pi in p]}") + print(f" Zeros: {z}") + print(f" Controllable: {ctrl}") + print(f" Observable: {obs}") + + # --- Frequency response --- + omega, mag_db, phase_deg = frequency_response(ls.A, ls.B, ls.C, ls.D) + print(f"\n Bode data: {len(omega)} frequency points") + + # --- PID simulation --- + print("\n--- PID step response (Kp=10, Ki=5, Kd=0.5) ---") + times, theta = simulate_pid(kp=10.0, ki=5.0, kd=0.5) + metrics = step_response_metrics(times, theta, setpoint=1.0) + print(f" Rise time: {metrics.rise_time:.3f} s") + print(f" Settling: {metrics.settling_time:.3f} s") + print(f" Overshoot: {metrics.overshoot_pct:.1f}%") + print(f" SS error: {metrics.steady_state_error:.4f}") + print(" (Integral action eliminates steady-state error)") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson4_disturbance.py b/packages/gds-examples/gds_examples/control_tutorial/lesson4_disturbance.py new file mode 100644 index 0000000..ff665d7 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson4_disturbance.py @@ -0,0 +1,160 @@ +"""Lesson 4 -- Disturbance Rejection and Sensitivity Analysis. + +Add a load torque disturbance to the motor and analyze how the +feedback loop handles it. Introduce feedforward compensation and +the Gang of Six sensitivity functions that quantify the trade-offs +between reference tracking, disturbance rejection, and noise immunity. + +New concepts: + - 2-input plant (voltage V + load torque tau_load) + - Feedforward disturbance compensation + - sensitivity() -> Gang of Six (S, T, CS, PS, KS, KPS) + - Sensitivity function Bode plots + +GDS Decomposition: + X = (theta, omega) + U = (theta_ref, tau_load) + g = pid_controller + feedforward + f = motor_dynamics (with disturbance) + Theta = {J, b, K_t, Kp, Ki, Kd} + +Composition: + (reference | disturbance | sensor) >> (controller | feedforward) + >> motor_dynamics .loop(dynamics -> sensor) +""" + +from gds_domains.symbolic.transfer import TransferFunction, sensitivity + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Transfer function helpers +# --------------------------------------------------------------------------- + + +def get_plant_tf() -> TransferFunction: + """Open-loop plant: P(s) = K_t/J / (s^2 + (b/J)*s). + + H(s) = (K_t/J) / (s^2 + (b/J)*s) = 1.0 / (s^2 + 10*s) + """ + return TransferFunction( + num=[K_t / J], + den=[1.0, b / J, 0.0], # s^2 + 10s + ) + + +def get_pid_tf(kp: float = 10.0, ki: float = 5.0, kd: float = 0.5) -> TransferFunction: + """PID controller: K(s) = Kp + Ki/s + Kd*s = (Kd*s^2 + Kp*s + Ki) / s.""" + return TransferFunction( + num=[kd, kp, ki], # Kd*s^2 + Kp*s + Ki + den=[1.0, 0.0], # s + ) + + +def analyze_sensitivity( + kp: float = 10.0, ki: float = 5.0, kd: float = 0.5 +) -> dict[str, TransferFunction]: + """Compute the Gang of Six sensitivity functions. + + Returns dict with keys S, T, CS, PS, KS, KPS. + """ + plant = get_plant_tf() + controller = get_pid_tf(kp, ki, kd) + return sensitivity(plant, controller) + + +def simulate_disturbance( + kp: float = 10.0, + ki: float = 5.0, + kd: float = 0.5, + tau_load: float = 0.005, + use_feedforward: bool = False, + t_end: float = 10.0, +) -> tuple[list[float], list[float]]: + """Simulate motor with step disturbance at t=t_end/2. + + The disturbance tau_load is applied as a constant after halfway. + If use_feedforward=True, adds V_ff = tau_load / K_t to cancel it. + + Returns (times, theta_values). + """ + from gds_continuous import ODEModel, ODESimulation + + t_disturb = t_end / 2 + + def motor_with_disturbance(t, y, params): + theta, omega, err_int = y + ref = params.get("theta_ref", 1.0) + kp_ = params.get("Kp", kp) + ki_ = params.get("Ki", ki) + kd_ = params.get("Kd", kd) + + error = ref - theta + v_fb = kp_ * error + ki_ * err_int + kd_ * (0 - omega) + + # Disturbance: step at t_disturb + d = tau_load if t >= t_disturb else 0.0 + + # Feedforward: cancel known disturbance + v_ff = (d / K_t) if use_feedforward else 0.0 + + v = v_fb + v_ff + + dtheta = omega + domega = -(b / J) * omega + (K_t / J) * v - (1 / J) * d + derr = error + return [dtheta, domega, derr] + + ode_model = ODEModel( + state_names=["theta", "omega", "error_int"], + initial_state={"theta": 0.0, "omega": 0.0, "error_int": 0.0}, + rhs=motor_with_disturbance, + params={"theta_ref": [1.0], "Kp": [kp], "Ki": [ki], "Kd": [kd]}, + ) + sim = ODESimulation(model=ode_model, t_span=(0.0, t_end), solver="RK45") + results = sim.run() + return results.times, results.state_array("theta") + + +if __name__ == "__main__": + print("=== Lesson 4: Disturbance Rejection + Gang of Six ===\n") + + # --- Gang of Six --- + gang = analyze_sensitivity() + print("Gang of Six sensitivity functions:") + for name, tf in gang.items(): + print( + f" {name:>4}: num={[f'{c:.3f}' for c in tf.num]}, " + f"den={[f'{c:.3f}' for c in tf.den]}" + ) + + # Evaluate S(0) and T(0) — DC sensitivity + s_dc = gang["S"].num[-1] / gang["S"].den[-1] + t_dc = gang["T"].num[-1] / gang["T"].den[-1] + print(f"\n S(0) = {s_dc:.4f} (ideal: 0 for perfect disturbance rejection)") + print(f" T(0) = {t_dc:.4f} (ideal: 1 for perfect reference tracking)") + print(f" S(0) + T(0) = {s_dc + t_dc:.4f} (must equal 1)") + + # --- Disturbance simulation --- + print("\n--- Step disturbance at t=5s (tau_load=0.005 N*m) ---") + t_fb, y_fb = simulate_disturbance(use_feedforward=False) + t_ff, y_ff = simulate_disturbance(use_feedforward=True) + + # Measure disturbance effect (deviation from setpoint after disturbance) + mid = len(t_fb) // 2 + max_dev_fb = max(abs(v - 1.0) for v in y_fb[mid:]) + max_dev_ff = max(abs(v - 1.0) for v in y_ff[mid:]) + + print(f" FB-only: max deviation = {max_dev_fb:.4f} rad") + print(f" FF+FB: max deviation = {max_dev_ff:.4f} rad") + print( + f" Feedforward reduces disturbance effect by " + f"{(1 - max_dev_ff / max_dev_fb) * 100:.0f}%" + ) diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson5_delay.py b/packages/gds-examples/gds_examples/control_tutorial/lesson5_delay.py new file mode 100644 index 0000000..9fc9861 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson5_delay.py @@ -0,0 +1,164 @@ +"""Lesson 5 -- Sensor Delay, Pade Approximation, and Stability Margins. + +Real sensors have sampling delay, communication latency, and processing +time. We model this as a pure time delay e^{-s*tau}, approximate it +with a Pade rational polynomial, and analyze how the delay degrades +gain and phase margins. The Nyquist plot shows how close the delayed +system comes to encircling the critical point (-1, 0). + +New concepts: + - pade_approximation() (from gds_domains.symbolic.delay) + - delay_system() for cascading delay with plant TF + - gain_margin(), phase_margin() (from gds_analysis.linear) + - nyquist_plot(), root_locus_plot() (from gds_viz.frequency) + - Stability margin degradation analysis + +GDS Decomposition: + X = (theta, omega) + U = theta_ref + g = pid_controller (with delayed measurement) + f = motor_dynamics + Theta = {J, b, K_t, Kp, Ki, Kd, tau_delay} + +Composition: + reference >> pid_controller >> motor_dynamics >> [DELAY] >> sensor + .loop(sensor -> pid_controller) +""" + +from gds_domains.symbolic.delay import pade_approximation +from gds_domains.symbolic.transfer import TransferFunction, _tf_multiply + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Transfer function construction +# --------------------------------------------------------------------------- + + +def get_plant_tf() -> TransferFunction: + """Open-loop plant: P(s) = (K_t/J) / (s^2 + (b/J)*s).""" + return TransferFunction(num=[K_t / J], den=[1.0, b / J, 0.0]) + + +def get_pid_tf(kp: float = 10.0, ki: float = 5.0, kd: float = 0.5) -> TransferFunction: + """PID controller TF: K(s) = (Kd*s^2 + Kp*s + Ki) / s.""" + return TransferFunction(num=[kd, kp, ki], den=[1.0, 0.0]) + + +def get_loop_tf(delay: float = 0.0, pade_order: int = 3) -> TransferFunction: + """Open-loop transfer function L(s) = K(s) * P(s) [* Pade(s)]. + + If delay > 0, cascades the Pade approximation of e^{-s*tau}. + """ + plant = get_plant_tf() + controller = get_pid_tf() + loop = _tf_multiply(controller, plant) + + if delay > 0: + pade_tf = pade_approximation(delay, order=pade_order) + loop = _tf_multiply(loop, pade_tf) + + return loop + + +def compare_margins( + delays: list[float] | None = None, +) -> list[dict]: + """Compute gain/phase margins for multiple delay values. + + Returns list of dicts with keys: delay, gm_db, gm_freq, pm_deg, pm_freq. + """ + from gds_analysis.linear import gain_margin, phase_margin + + if delays is None: + delays = [0.0, 0.02, 0.05, 0.1] + + results = [] + for tau in delays: + loop = get_loop_tf(delay=tau) + gm_db, gm_freq = gain_margin(loop.num, loop.den) + pm_deg, pm_freq = phase_margin(loop.num, loop.den) + results.append( + { + "delay": tau, + "gm_db": gm_db, + "gm_freq": gm_freq, + "pm_deg": pm_deg, + "pm_freq": pm_freq, + } + ) + return results + + +def get_nyquist_data( + delay: float = 0.05, +) -> tuple[list[float], list[float]]: + """Compute Nyquist plot data for the delayed loop TF. + + Returns (real_parts, imag_parts) of L(jw). + """ + + loop = get_loop_tf(delay=delay) + + # Build a state-space from the TF for frequency_response + # For simplicity, use the TF coefficients directly + import numpy as np + + w = np.logspace(-2, 3, 2000) + s_vals = 1j * w + + # Evaluate L(jw) directly from polynomial coefficients + num = np.array(loop.num) + den = np.array(loop.den) + + def eval_poly(coeffs, s): + n = len(coeffs) - 1 + return sum(c * s ** (n - i) for i, c in enumerate(coeffs)) + + H = np.array([eval_poly(num, s) / eval_poly(den, s) for s in s_vals]) + + return np.real(H).tolist(), np.imag(H).tolist() + + +if __name__ == "__main__": + print("=== Lesson 5: Sensor Delay and Stability Margins ===\n") + + # --- Pade approximation --- + pade = pade_approximation(0.05, order=3) + print("Pade(tau=0.05s, order=3):") + print(f" num = {[f'{c:.6f}' for c in pade.num]}") + print(f" den = {[f'{c:.6f}' for c in pade.den]}") + + # --- Margin comparison --- + print( + f"\n{'Delay':>8} {'GM (dB)':>10} {'GM freq':>10} " + f"{'PM (deg)':>10} {'PM freq':>10}" + ) + print("-" * 58) + + margins = compare_margins() + for m in margins: + gm = f"{m['gm_db']:.1f}" if m["gm_db"] != float("inf") else "inf" + pm = f"{m['pm_deg']:.1f}" if m["pm_deg"] != float("inf") else "inf" + gf = f"{m['gm_freq']:.2f}" if m["gm_freq"] == m["gm_freq"] else "n/a" + pf = f"{m['pm_freq']:.2f}" if m["pm_freq"] == m["pm_freq"] else "n/a" + print(f"{m['delay']:8.3f} {gm:>10} {gf:>10} {pm:>10} {pf:>10}") + + # --- Stability warning --- + for m in margins: + if m["pm_deg"] != float("inf") and m["pm_deg"] < 30: + print( + f"\n WARNING: Phase margin at delay={m['delay']}s is " + f"{m['pm_deg']:.1f} deg — below 30 deg safety threshold!" + ) + + print("\n As delay increases, phase margin shrinks.") + print(" At some point the system becomes unstable.") + print(" This is why sensor bandwidth matters in control design.") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson6_lqr.py b/packages/gds-examples/gds_examples/control_tutorial/lesson6_lqr.py new file mode 100644 index 0000000..5822540 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson6_lqr.py @@ -0,0 +1,199 @@ +"""Lesson 6 -- LQR Optimal Control vs PID. + +Replace the hand-tuned PID with a Linear Quadratic Regulator (LQR). +LQR computes the optimal gain matrix K by minimizing a cost function +J = integral(x'Qx + u'Ru) dt, balancing state error against control +effort. We compare step responses, design a Kalman observer, and +compute a gain schedule for varying load inertia. + +New concepts: + - lqr() -> (K, P, E) from gds_analysis.linear + - kalman() -> (L, P) for observer design + - gain_schedule() for multiple operating points + - Q/R cost matrix design philosophy + - compare_responses() for LQR vs PID + +GDS Decomposition: + X = (theta, omega) + U = theta_ref + g = lqr_controller: V = -K @ [theta - theta_ref, omega] + f = motor_dynamics + Theta = {J, b, K_t, Q, R} + +Composition: + (reference | sensor) >> lqr_controller >> motor_dynamics + .loop(dynamics -> sensor) +""" + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Plant linearization (analytical, matching Lesson 3) +# --------------------------------------------------------------------------- + + +def get_plant_ss() -> tuple[ + list[list[float]], + list[list[float]], + list[list[float]], + list[list[float]], +]: + """Return (A, B, C, D) for the open-loop DC motor plant.""" + A = [[0.0, 1.0], [0.0, -b / J]] + B = [[0.0], [K_t / J]] + C = [[1.0, 0.0]] + D = [[0.0]] + return A, B, C, D + + +# --------------------------------------------------------------------------- +# LQR design +# --------------------------------------------------------------------------- + + +def design_lqr( + q_diag: tuple[float, float] = (100.0, 1.0), + r_val: float = 1.0, +) -> tuple[list[list[float]], list[list[float]], list[complex]]: + """Design LQR gains for the DC motor. + + Args: + q_diag: Diagonal of Q (position error weight, velocity error weight). + r_val: Control effort weight R. + + Returns (K, P, E): gain, Riccati solution, closed-loop eigenvalues. + """ + from gds_analysis.linear import lqr + + A, B, _, _ = get_plant_ss() + Q = [[q_diag[0], 0.0], [0.0, q_diag[1]]] + R = [[r_val]] + return lqr(A, B, Q, R) + + +def design_kalman() -> tuple[list[list[float]], list[list[float]]]: + """Design Kalman observer for the DC motor. + + Returns (L, P): observer gain, error covariance. + """ + from gds_analysis.linear import kalman + + A, _, C, _ = get_plant_ss() + Q_process = [[0.01, 0.0], [0.0, 0.1]] # process noise covariance + R_measurement = [[0.001]] # measurement noise covariance + return kalman(A, C, Q_process, R_measurement) + + +def simulate_lqr( + K: list[list[float]], + theta_ref: float = 1.0, + t_end: float = 3.0, +) -> tuple[list[float], list[float]]: + """Simulate LQR-controlled motor: V = -K @ [theta - ref, omega]. + + Returns (times, theta_values). + """ + from gds_continuous import ODEModel, ODESimulation + + k1, k2 = K[0][0], K[0][1] + + def lqr_ode(t, y, params): + theta, omega = y + ref = params.get("theta_ref", 1.0) + v = -k1 * (theta - ref) - k2 * omega + dtheta = omega + domega = -(b / J) * omega + (K_t / J) * v + return [dtheta, domega] + + ode_model = ODEModel( + state_names=["theta", "omega"], + initial_state={"theta": 0.0, "omega": 0.0}, + rhs=lqr_ode, + params={"theta_ref": [theta_ref]}, + ) + sim = ODESimulation(model=ode_model, t_span=(0.0, t_end), solver="RK45") + results = sim.run() + return results.times, results.state_array("theta") + + +def run_gain_schedule() -> list[tuple[dict, list[list[float]], list[complex]]]: + """Compute LQR gains at 3 different load inertias. + + Simulates the motor with different payloads attached. + """ + from gds_analysis.linear import gain_schedule + + def linearize_at_inertia(point): + j = point["J"] + A = [[0.0, 1.0], [0.0, -b / j]] + B = [[0.0], [K_t / j]] + C = [[1.0, 0.0]] + D = [[0.0]] + return A, B, C, D + + points = [{"J": 0.005}, {"J": 0.01}, {"J": 0.02}] + Q = [[100.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + + return gain_schedule(linearize_at_inertia, points, Q, R) + + +if __name__ == "__main__": + from gds_analysis.response import step_response_metrics + + print("=== Lesson 6: LQR Optimal Control ===\n") + + # --- LQR design --- + K, P, E = design_lqr() + print(f"LQR gain K = [{K[0][0]:.4f}, {K[0][1]:.4f}]") + print(f"Riccati P diagonal = [{P[0][0]:.4f}, {P[1][1]:.4f}]") + print("Closed-loop eigenvalues:") + for e in E: + print(f" {e.real:.4f} + {e.imag:.4f}j") + print(f"Closed-loop stable: {all(e.real < 0 for e in E)}") + + # --- LQR step response --- + print("\n--- LQR step response ---") + t_lqr, y_lqr = simulate_lqr(K) + m_lqr = step_response_metrics(t_lqr, y_lqr, setpoint=1.0) + print(f" Rise time: {m_lqr.rise_time:.3f} s") + print(f" Settling: {m_lqr.settling_time:.3f} s") + print(f" Overshoot: {m_lqr.overshoot_pct:.1f}%") + + # --- Compare with PID (from Lesson 3) --- + from gds_examples.control_tutorial.lesson3_pid_transfer import simulate_pid + + t_pid, y_pid = simulate_pid(kp=10.0, ki=5.0, kd=0.5, t_end=3.0) + m_pid = step_response_metrics(t_pid, y_pid, setpoint=1.0) + + print("\n--- LQR vs PID comparison ---") + print(f"{'Metric':>15} {'LQR':>8} {'PID':>8}") + print("-" * 36) + print(f"{'Rise time':>15} {m_lqr.rise_time:8.3f} {m_pid.rise_time:8.3f}") + print(f"{'Settling':>15} {m_lqr.settling_time:8.3f} {m_pid.settling_time:8.3f}") + os_lqr, os_pid = m_lqr.overshoot_pct, m_pid.overshoot_pct + print(f"{'Overshoot %':>15} {os_lqr:8.1f} {os_pid:8.1f}") + + # --- Kalman observer --- + print("\n--- Kalman observer ---") + L, P_obs = design_kalman() + print(f" Observer gain L = [{L[0][0]:.4f}, {L[1][0]:.4f}]'") + + # --- Gain schedule --- + print("\n--- Gain schedule (varying inertia) ---") + schedule = run_gain_schedule() + print(f"{'J':>8} {'K1':>10} {'K2':>10} {'eig1':>12} {'eig2':>12}") + print("-" * 58) + for point, k, e in schedule: + print( + f"{point['J']:8.3f} {k[0][0]:10.4f} {k[0][1]:10.4f}" + f" {e[0].real:12.4f} {e[1].real:12.4f}" + ) + print("\n Gains change with inertia — one K doesn't fit all loads.") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson7_discrete.py b/packages/gds-examples/gds_examples/control_tutorial/lesson7_discrete.py new file mode 100644 index 0000000..3d9d647 --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson7_discrete.py @@ -0,0 +1,220 @@ +"""Lesson 7 -- Discretization and Discrete-Time Simulation. + +Continuous-time controllers run on digital microcontrollers that update +at fixed sample rates. We convert the plant and LQR controller from +continuous to discrete time using Tustin and ZOH methods, verify +discrete stability (eigenvalues inside the unit circle), and simulate +in gds-sim's discrete-time engine. + +New concepts: + - discretize() with Tustin / ZOH methods (from gds_analysis.linear) + - is_stable(A, continuous=False) for discrete stability + - dlqr() for discrete-time LQR + - gds_sim Model / Simulation for timestep-based execution + - Sample rate sweep: effect of dt on performance + +GDS Decomposition: + X[k] = (theta[k], omega[k]) + U[k] = theta_ref + g[k] = discrete_lqr: V[k] = -K_d @ (X[k] - X_ref) + f[k] = X[k+1] = Ad @ X[k] + Bd @ V[k] + Theta = {Ad, Bd, K_d, dt} + +Composition: + reference >> discrete_controller >> discrete_dynamics + .loop(dynamics -> sensor) [temporal, COVARIANT] +""" + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Discretization +# --------------------------------------------------------------------------- + + +def discretize_plant( + dt: float = 0.01, + method: str = "tustin", +) -> tuple[ + list[list[float]], + list[list[float]], + list[list[float]], + list[list[float]], +]: + """Discretize the DC motor plant at sample period dt. + + Returns (Ad, Bd, Cd, Dd). + """ + from gds_analysis.linear import discretize + + A = [[0.0, 1.0], [0.0, -b / J]] + B = [[0.0], [K_t / J]] + C = [[1.0, 0.0]] + D = [[0.0]] + return discretize(A, B, C, D, dt, method=method) + + +def design_dlqr( + dt: float = 0.01, + q_diag: tuple[float, float] = (100.0, 1.0), + r_val: float = 1.0, +) -> tuple[list[list[float]], list[list[float]], list[complex]]: + """Design discrete-time LQR for the discretized plant. + + Returns (Kd, Pd, Ed): discrete gain, Riccati solution, eigenvalues. + """ + from gds_analysis.linear import dlqr + + Ad, Bd, _Cd, _Dd = discretize_plant(dt, method="zoh") + Q = [[q_diag[0], 0.0], [0.0, q_diag[1]]] + R = [[r_val]] + return dlqr(Ad, Bd, Q, R) + + +# --------------------------------------------------------------------------- +# Discrete simulation via gds-sim +# --------------------------------------------------------------------------- + + +def simulate_discrete( + dt: float = 0.01, + theta_ref: float = 1.0, + n_steps: int = 300, +) -> tuple[list[float], list[float]]: + """Simulate discrete-time LQR control using gds-sim. + + Returns (times, theta_values). + """ + from gds_sim import Model, Simulation + + Ad, Bd, _Cd, _Dd = discretize_plant(dt, method="zoh") + Kd, _Pd, _Ed = design_dlqr(dt) + + k1, k2 = Kd[0][0], Kd[0][1] + + def controller_policy(_params, _substep, _history, state): + theta = state["theta"] + omega = state["omega"] + v = -k1 * (theta - theta_ref) - k2 * omega + return {"voltage": v} + + def theta_suf(_params, _substep, _history, state, signal): + theta = state["theta"] + omega = state["omega"] + v = signal["voltage"] + theta_next = Ad[0][0] * theta + Ad[0][1] * omega + Bd[0][0] * v + return ("theta", theta_next) + + def omega_suf(_params, _substep, _history, state, signal): + theta = state["theta"] + omega = state["omega"] + v = signal["voltage"] + omega_next = Ad[1][0] * theta + Ad[1][1] * omega + Bd[1][0] * v + return ("omega", omega_next) + + model = Model( + initial_state={"theta": 0.0, "omega": 0.0}, + state_update_blocks=[ + { + "policies": {"controller": controller_policy}, + "variables": {"theta": theta_suf, "omega": omega_suf}, + }, + ], + params={}, + ) + + sim = Simulation(model=model, timesteps=n_steps, runs=1) + results = sim.run() + rows = results.to_list() + + times = [i * dt for i in range(len(rows))] + theta_vals = [row["theta"] for row in rows] + + return times, theta_vals + + +# --------------------------------------------------------------------------- +# Sample rate sweep +# --------------------------------------------------------------------------- + + +def sample_rate_sweep( + rates: list[float] | None = None, +) -> list[dict]: + """Sweep sample rates and collect settling time metrics. + + Returns list of dicts with keys: dt, settling_time, overshoot_pct, stable. + """ + from gds_analysis.linear import is_stable + from gds_analysis.response import step_response_metrics + + if rates is None: + rates = [0.001, 0.005, 0.01, 0.02, 0.05] + + results = [] + for dt in rates: + Ad, _Bd, _Cd, _Dd = discretize_plant(dt, method="zoh") + stable = is_stable(Ad, continuous=False) + + times, theta = simulate_discrete(dt=dt, n_steps=int(3.0 / dt)) + metrics = step_response_metrics(times, theta, setpoint=1.0) + + results.append( + { + "dt": dt, + "settling_time": metrics.settling_time, + "overshoot_pct": metrics.overshoot_pct, + "stable": stable, + } + ) + return results + + +if __name__ == "__main__": + from gds_analysis.linear import eigenvalues + + print("=== Lesson 7: Discretization + gds-sim ===\n") + + # --- Discretize --- + dt = 0.01 + Ad, Bd, Cd, Dd = discretize_plant(dt, method="zoh") + print(f"Discrete plant (ZOH, dt={dt}s):") + print(f" Ad = {Ad}") + print(f" Bd = {Bd}") + + eigs = eigenvalues(Ad) + print(f" Eigenvalues: {[f'{abs(e):.6f}' for e in eigs]}") + print(f" |lambda| < 1: {all(abs(e) < 1.0 for e in eigs)}") + + # --- Discrete LQR --- + Kd, _Pd, Ed = design_dlqr(dt) + print(f"\nDiscrete LQR gain Kd = [{Kd[0][0]:.4f}, {Kd[0][1]:.4f}]") + print(f" Closed-loop |lambda|: {[f'{abs(e):.6f}' for e in Ed]}") + + # --- Discrete simulation --- + print(f"\n--- Discrete simulation (dt={dt}s, 300 steps) ---") + times, theta = simulate_discrete(dt=dt) + from gds_analysis.response import step_response_metrics + + m = step_response_metrics(times, theta, setpoint=1.0) + print(f" Rise time: {m.rise_time:.3f} s") + print(f" Settling: {m.settling_time:.3f} s") + print(f" Overshoot: {m.overshoot_pct:.1f}%") + + # --- Sample rate sweep --- + print("\n--- Sample rate sweep ---") + print(f"{'dt':>8} {'Settle':>8} {'OS%':>8} {'Stable':>8}") + print("-" * 38) + for r in sample_rate_sweep(): + print( + f"{r['dt']:8.3f} {r['settling_time']:8.3f} " + f"{r['overshoot_pct']:8.1f} {'yes' if r['stable'] else 'NO':>8}" + ) + print("\n Slower sample rate -> worse performance, eventually unstable.") diff --git a/packages/gds-examples/gds_examples/control_tutorial/lesson8_lyapunov.py b/packages/gds-examples/gds_examples/control_tutorial/lesson8_lyapunov.py new file mode 100644 index 0000000..88fc92d --- /dev/null +++ b/packages/gds-examples/gds_examples/control_tutorial/lesson8_lyapunov.py @@ -0,0 +1,182 @@ +"""Lesson 8 -- Lyapunov Stability Proof and Passivity Certificate. + +Move from numerical stability checks (eigenvalues) to formal symbolic +proofs. A Lyapunov function V(x) > 0 with dV/dt < 0 is the control +theorist's invariant — it guarantees that energy is always dissipating. +We find V for the closed-loop motor, verify it symbolically, and then +prove passivity: the motor dissipates more energy than it stores. + +New concepts: + - quadratic_lyapunov(P, A) -- verify V = x'Px via eigenvalue check + - find_quadratic_lyapunov(A) -- solve Lyapunov equation A'P + PA = -Q + - lyapunov_candidate(V, f) -- verify custom V(x) symbolically + - passivity_certificate(V, s, f, h) -- energy dissipation proof + - LyapunovResult / PassivityResult result types + +GDS Decomposition: + X = (theta, omega) + U = theta_ref + g = lqr_controller + f = motor_dynamics (closed-loop: A_cl = A - B*K) + Theta = {J, b, K_t, K_lqr} + +The Lyapunov function V(x) = x'Px is the "generalized energy" of the +closed-loop system. If V always decreases, the system is stable. +""" + +# --------------------------------------------------------------------------- +# Motor parameters +# --------------------------------------------------------------------------- + +J = 0.01 +b = 0.1 +K_t = 0.01 + + +# --------------------------------------------------------------------------- +# Closed-loop system +# --------------------------------------------------------------------------- + + +def get_closed_loop_A() -> list[list[float]]: + """Compute closed-loop A_cl = A - B*K using LQR from Lesson 6. + + Returns the 2x2 closed-loop state matrix. + """ + from gds_analysis.linear import lqr + + A = [[0.0, 1.0], [0.0, -b / J]] + B = [[0.0], [K_t / J]] + Q = [[100.0, 0.0], [0.0, 1.0]] + R = [[1.0]] + + K, _P, _E = lqr(A, B, Q, R) + + # A_cl = A - B*K + A_cl = [ + [A[0][0] - B[0][0] * K[0][0], A[0][1] - B[0][0] * K[0][1]], + [A[1][0] - B[1][0] * K[0][0], A[1][1] - B[1][0] * K[0][1]], + ] + return A_cl + + +# --------------------------------------------------------------------------- +# Lyapunov analysis +# --------------------------------------------------------------------------- + + +def prove_quadratic_lyapunov(): + """Find and verify a quadratic Lyapunov function for the closed-loop motor. + + Returns (P, result) where P is the Lyapunov matrix and result + is a LyapunovResult with status fields. + """ + from gds_proof.analysis.lyapunov import find_quadratic_lyapunov + + A_cl = get_closed_loop_A() + return find_quadratic_lyapunov(A_cl) + + +def verify_custom_lyapunov(): + """Verify V(x) = theta^2 + omega^2 as a Lyapunov candidate. + + This is a simple "energy-like" function (not the optimal one). + It may or may not prove stable depending on the dynamics. + + Returns LyapunovResult. + """ + from gds_proof.analysis.lyapunov import lyapunov_candidate + + A_cl = get_closed_loop_A() + + # Build state transition from closed-loop A + # d(theta)/dt = A_cl[0][0]*theta + A_cl[0][1]*omega + # d(omega)/dt = A_cl[1][0]*theta + A_cl[1][1]*omega + return lyapunov_candidate( + V_expr="theta**2 + omega**2", + state_transition={ + "theta": f"{A_cl[0][0]}*theta + {A_cl[0][1]}*omega", + "omega": f"{A_cl[1][0]}*theta + {A_cl[1][1]}*omega", + }, + state_symbols=["theta", "omega"], + continuous=True, + ) + + +def prove_passivity(): + """Prove passivity of the open-loop motor plant. + + Storage function: V = 0.5 * J * omega^2 (kinetic energy) + Supply rate: s = K_t * V_input * omega (torque times velocity = power in) + Dynamics: d(omega)/dt = -(b/J)*omega + (K_t/J)*V_input + + dV/dt = J * omega * domega/dt = -b*omega^2 + K_t*V_input*omega + Dissipation: dV/dt - s = -b*omega^2 <= 0 + The damping term b*omega^2 is always non-negative, so the motor + is passive (it dissipates energy through friction). + + Returns PassivityResult. + """ + from gds_proof.analysis.lyapunov import passivity_certificate + + return passivity_certificate( + V_expr=f"{J}/2 * omega**2", + supply_rate=f"{K_t} * V_input * omega", + state_transition={ + "omega": f"-({b}/{J})*omega + ({K_t}/{J})*V_input", + }, + output_map={"y_omega": "omega"}, + state_symbols=["omega"], + input_symbols=["V_input"], + ) + + +if __name__ == "__main__": + from gds_analysis.linear import eigenvalues + + print("=== Lesson 8: Lyapunov Stability Proof ===\n") + + # --- Closed-loop system --- + A_cl = get_closed_loop_A() + eigs = eigenvalues(A_cl) + print("Closed-loop A_cl:") + print(f" [{A_cl[0][0]:8.4f}, {A_cl[0][1]:8.4f}]") + print(f" [{A_cl[1][0]:8.4f}, {A_cl[1][1]:8.4f}]") + print(f" Eigenvalues: {[f'{e.real:.4f}' for e in eigs]}") + + # --- Find quadratic Lyapunov --- + print("\n--- Quadratic Lyapunov (find P s.t. A'P + PA = -I) ---") + result = prove_quadratic_lyapunov() + if result is not None: + P, lyap = result + print(f" P = [[{P[0][0]:.4f}, {P[0][1]:.4f}],") + print(f" [{P[1][0]:.4f}, {P[1][1]:.4f}]]") + print(f" Positive definite: {lyap.positive_definite}") + print(f" Decreasing: {lyap.decreasing}") + print(f" STABLE: {lyap.stable}") + else: + print(" Failed to find P (system may be unstable)") + + # --- Custom Lyapunov candidate --- + print("\n--- Custom candidate: V = theta^2 + omega^2 ---") + custom = verify_custom_lyapunov() + print(f" Positive definite: {custom.positive_definite}") + print(f" Decreasing: {custom.decreasing}") + print(f" Stable: {custom.stable}") + if custom.dV_expr is not None: + print(f" dV/dt = {custom.dV_expr}") + + # --- Passivity --- + print("\n--- Passivity certificate (open-loop motor) ---") + passivity = prove_passivity() + print(f" Storage: V = {passivity.storage_function}") + print(f" Supply rate: s = {passivity.supply_rate}") + print(f" Dissipation: {passivity.dissipation_proved}") + print(f" Passive: {passivity.passive}") + print("\n Proof trace:") + for line in passivity.details: + print(f" {line}") + + print("\n The motor is passive: friction dissipates energy.") + print(" Any passive controller guarantees closed-loop stability,") + print(" even with unmodeled dynamics (passivity theorem).") diff --git a/packages/gds-examples/tests/test_control_tutorial.py b/packages/gds-examples/tests/test_control_tutorial.py new file mode 100644 index 0000000..77e1603 --- /dev/null +++ b/packages/gds-examples/tests/test_control_tutorial.py @@ -0,0 +1,331 @@ +"""Tests for the DC Motor Control Tutorial (Lessons 1-8).""" + +import pytest + +# --------------------------------------------------------------------------- +# Lesson 1: Plant model +# --------------------------------------------------------------------------- + + +class TestLesson1Plant: + def test_model_has_two_states(self) -> None: + from gds_examples.control_tutorial.lesson1_plant import build_model + + model = build_model() + assert len(model.states) == 2 + assert {s.name for s in model.states} == {"theta", "omega"} + + def test_model_has_one_input(self) -> None: + from gds_examples.control_tutorial.lesson1_plant import build_model + + model = build_model() + assert len(model.inputs) == 1 + assert model.inputs[0].name == "V" + + def test_spec_validates(self) -> None: + from gds_examples.control_tutorial.lesson1_plant import build_spec + + spec = build_spec() + errors = spec.validate_spec() + assert errors == [] + + def test_system_compiles(self) -> None: + from gds_examples.control_tutorial.lesson1_plant import build_system + + system = build_system() + assert len(system.blocks) > 0 + + def test_open_loop_simulation_runs(self) -> None: + from gds_examples.control_tutorial.lesson1_plant import simulate_open_loop + + times, theta, omega = simulate_open_loop(voltage=1.0, t_end=1.0) + assert len(times) > 10 + # Motor should spin up: omega > 0 at end + assert omega[-1] > 0 + # Position should increase + assert theta[-1] > 0 + + +# --------------------------------------------------------------------------- +# Lesson 2: P control +# --------------------------------------------------------------------------- + + +class TestLesson2PControl: + def test_model_builds(self) -> None: + from gds_examples.control_tutorial.lesson2_p_control import build_model + + model = build_model() + assert len(model.states) == 2 + assert "Kp" in model.symbolic_params + + def test_step_response_converges(self) -> None: + from gds_examples.control_tutorial.lesson2_p_control import ( + simulate_p_control, + ) + + _times, theta = simulate_p_control(kp=5.0, theta_ref=1.0, t_end=5.0) + # Should approach setpoint (with some SSE for P control) + assert theta[-1] > 0.5 + + def test_step_metrics_reasonable(self) -> None: + from gds_examples.control_tutorial.lesson2_p_control import ( + analyze_step_response, + ) + + _t, _y, metrics = analyze_step_response(kp=5.0) + assert metrics.rise_time > 0 + assert metrics.settling_time > 0 + assert metrics.overshoot_pct >= 0 + + +# --------------------------------------------------------------------------- +# Lesson 3: PID + transfer functions +# --------------------------------------------------------------------------- + + +class TestLesson3PIDTransfer: + def test_linearization_dimensions(self) -> None: + from gds_examples.control_tutorial.lesson3_pid_transfer import ( + get_plant_linearization, + ) + + ls = get_plant_linearization() + assert len(ls.A) == 2 + assert len(ls.A[0]) == 2 + assert len(ls.B) == 2 + assert len(ls.B[0]) == 1 + + def test_transfer_function_order(self) -> None: + from gds_examples.control_tutorial.lesson3_pid_transfer import ( + analyze_transfer_function, + ) + + tf, _p, _z, _ctrl, _obs = analyze_transfer_function() + # Plant TF denominator should be order 2 + assert len(tf.den) == 3 + + def test_plant_is_controllable(self) -> None: + from gds_examples.control_tutorial.lesson3_pid_transfer import ( + analyze_transfer_function, + ) + + _tf, _p, _z, ctrl, obs = analyze_transfer_function() + assert ctrl is True + assert obs is True + + def test_open_loop_poles_stable(self) -> None: + """Open-loop motor has poles at 0 and -b/J — marginally stable.""" + from gds_examples.control_tutorial.lesson3_pid_transfer import ( + analyze_transfer_function, + ) + + _tf, plant_poles, _z, _ctrl, _obs = analyze_transfer_function() + assert len(plant_poles) == 2 + # One pole at 0, one at -10 + reals = sorted([p.real for p in plant_poles]) + assert reals[0] == pytest.approx(-10.0, abs=0.5) + assert reals[1] == pytest.approx(0.0, abs=0.1) + + def test_pid_eliminates_sse(self) -> None: + from gds_analysis.response import step_response_metrics + from gds_examples.control_tutorial.lesson3_pid_transfer import simulate_pid + + times, theta = simulate_pid(kp=10.0, ki=5.0, kd=0.5, t_end=10.0) + metrics = step_response_metrics(times, theta, setpoint=1.0) + assert metrics.steady_state_error < 0.05 + + +# --------------------------------------------------------------------------- +# Lesson 4: Disturbance + sensitivity +# --------------------------------------------------------------------------- + + +class TestLesson4Disturbance: + def test_sensitivity_functions_exist(self) -> None: + from gds_examples.control_tutorial.lesson4_disturbance import ( + analyze_sensitivity, + ) + + gang = analyze_sensitivity() + assert set(gang.keys()) == {"S", "T", "CS", "PS", "KS", "KPS"} + + def test_s_plus_t_dc_equals_one(self) -> None: + from gds_examples.control_tutorial.lesson4_disturbance import ( + analyze_sensitivity, + ) + + gang = analyze_sensitivity() + s_dc = gang["S"].num[-1] / gang["S"].den[-1] + t_dc = gang["T"].num[-1] / gang["T"].den[-1] + assert s_dc + t_dc == pytest.approx(1.0, abs=1e-8) + + def test_feedforward_reduces_disturbance(self) -> None: + from gds_examples.control_tutorial.lesson4_disturbance import ( + simulate_disturbance, + ) + + _t_fb, y_fb = simulate_disturbance(use_feedforward=False) + _t_ff, y_ff = simulate_disturbance(use_feedforward=True) + + mid = len(y_fb) // 2 + dev_fb = max(abs(v - 1.0) for v in y_fb[mid:]) + dev_ff = max(abs(v - 1.0) for v in y_ff[mid:]) + assert dev_ff < dev_fb + + +# --------------------------------------------------------------------------- +# Lesson 5: Delay + margins +# --------------------------------------------------------------------------- + + +class TestLesson5Delay: + def test_pade_approximation_order(self) -> None: + from gds_domains.symbolic.delay import pade_approximation + + pade = pade_approximation(0.05, order=3) + assert len(pade.num) == 4 # order 3 -> degree 3 -> 4 coefficients + assert len(pade.den) == 4 + + def test_delay_reduces_margins(self) -> None: + from gds_examples.control_tutorial.lesson5_delay import compare_margins + + margins = compare_margins(delays=[0.0, 0.05]) + # No-delay margin should be better than delayed margin + no_delay = margins[0] + with_delay = margins[1] + + # At least one margin should degrade + if no_delay["pm_deg"] != float("inf") and with_delay["pm_deg"] != float("inf"): + assert with_delay["pm_deg"] < no_delay["pm_deg"] + + def test_nyquist_data_computable(self) -> None: + from gds_examples.control_tutorial.lesson5_delay import get_nyquist_data + + real, imag = get_nyquist_data(delay=0.05) + assert len(real) > 0 + assert len(imag) > 0 + + +# --------------------------------------------------------------------------- +# Lesson 6: LQR +# --------------------------------------------------------------------------- + + +class TestLesson6LQR: + def test_lqr_gain_shape(self) -> None: + from gds_examples.control_tutorial.lesson6_lqr import design_lqr + + K, _P, _E = design_lqr() + assert len(K) == 1 # single input + assert len(K[0]) == 2 # two states + + def test_closed_loop_stable(self) -> None: + from gds_examples.control_tutorial.lesson6_lqr import design_lqr + + _K, _P, E = design_lqr() + assert all(e.real < 0 for e in E) + + def test_lqr_step_response(self) -> None: + from gds_analysis.response import step_response_metrics + from gds_examples.control_tutorial.lesson6_lqr import ( + design_lqr, + simulate_lqr, + ) + + K, _P, _E = design_lqr() + times, theta = simulate_lqr(K, t_end=3.0) + metrics = step_response_metrics(times, theta, setpoint=1.0) + assert metrics.overshoot_pct < 50 + + def test_gain_schedule_varies(self) -> None: + from gds_examples.control_tutorial.lesson6_lqr import run_gain_schedule + + schedule = run_gain_schedule() + assert len(schedule) == 3 + # Velocity gains (K2) should differ across operating points + # (K1 position gain converges but K2 varies with inertia) + k2_values = [s[1][0][1] for s in schedule] + assert k2_values[0] != pytest.approx(k2_values[2], abs=0.01) + + +# --------------------------------------------------------------------------- +# Lesson 7: Discretization +# --------------------------------------------------------------------------- + + +class TestLesson7Discrete: + def test_dlqr_closed_loop_stable(self) -> None: + """Closed-loop discrete system should have eigenvalues inside unit circle.""" + from gds_examples.control_tutorial.lesson7_discrete import design_dlqr + + _Kd, _Pd, Ed = design_dlqr(dt=0.01) + assert all(abs(e) < 1.0 for e in Ed) + + def test_discrete_eigenvalues_inside_unit_circle(self) -> None: + from gds_examples.control_tutorial.lesson7_discrete import design_dlqr + + _Kd, _Pd, Ed = design_dlqr(dt=0.01) + assert all(abs(e) < 1.0 for e in Ed) + + def test_discrete_sim_runs(self) -> None: + from gds_examples.control_tutorial.lesson7_discrete import simulate_discrete + + times, theta = simulate_discrete(dt=0.01, n_steps=100) + assert len(times) > 50 + # Should approach setpoint + assert theta[-1] > 0.5 + + def test_dlqr_gain_computed(self) -> None: + from gds_examples.control_tutorial.lesson7_discrete import design_dlqr + + Kd, _Pd, Ed = design_dlqr(dt=0.01) + assert len(Kd) == 1 + assert len(Kd[0]) == 2 + assert all(abs(e) < 1.0 for e in Ed) + + +# --------------------------------------------------------------------------- +# Lesson 8: Lyapunov +# --------------------------------------------------------------------------- + + +class TestLesson8Lyapunov: + def test_find_lyapunov_succeeds(self) -> None: + from gds_examples.control_tutorial.lesson8_lyapunov import ( + prove_quadratic_lyapunov, + ) + + result = prove_quadratic_lyapunov() + assert result is not None + _P, lyap = result + assert lyap.stable is True + + def test_quadratic_lyapunov_proved(self) -> None: + from gds_examples.control_tutorial.lesson8_lyapunov import ( + prove_quadratic_lyapunov, + ) + + result = prove_quadratic_lyapunov() + assert result is not None + _P, lyap = result + assert lyap.positive_definite == "PROVED" + assert lyap.decreasing == "PROVED" + + def test_custom_candidate(self) -> None: + from gds_examples.control_tutorial.lesson8_lyapunov import ( + verify_custom_lyapunov, + ) + + result = verify_custom_lyapunov() + # V = theta^2 + omega^2 may or may not prove stable + # but should at least be positive definite + assert result.positive_definite == "PROVED" + assert result.dV_expr is not None + + def test_passivity_certificate(self) -> None: + from gds_examples.control_tutorial.lesson8_lyapunov import prove_passivity + + result = prove_passivity() + assert result.passive is True + assert result.dissipation_proved == "PROVED"