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/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/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/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..78ed64d --- /dev/null +++ b/packages/gds-analysis/tests/test_linear.py @@ -0,0 +1,340 @@ +"""Tests for numerical linear systems analysis.""" + +import pytest + +from gds_analysis.linear import ( + discretize, + dlqr, + eigenvalues, + frequency_response, + gain_margin, + gain_schedule, + 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 + + 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_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) + 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 +# --------------------------------------------------------------------------- + + +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_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", + ) + + +# --------------------------------------------------------------------------- +# 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 + + 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( + 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) + + +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-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..f4278c7 --- /dev/null +++ b/packages/gds-domains/tests/test_transfer.py @@ -0,0 +1,407 @@ +"""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, strict=False): + 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, 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: + 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"} + + 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-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" 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 new file mode 100644 index 0000000..04afe4c --- /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 + +Uses SymPy's simplification engine directly for positive-definiteness +and decrease condition checks via Hessian eigenvalue analysis. +""" + +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-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/__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..0250499 --- /dev/null +++ b/packages/gds-viz/gds_viz/frequency.py @@ -0,0 +1,452 @@ +"""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, + 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). + 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..ddfb750 --- /dev/null +++ b/packages/gds-viz/gds_viz/response.py @@ -0,0 +1,276 @@ +"""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 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] + 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" 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