gds-proof: Lyapunov stability proofs and passivity certificates

[rororowyourboat]

Summary

Add lyapunov.py to gds_proof/analysis/ providing Lyapunov stability proof templates, quadratic Lyapunov function verification, and passivity certificates. Extends the existing inductive safety proof engine with control-theoretic stability templates. SymPy only — no new dependencies.

Parent issue: #198 (classical control theory stack)

Motivation

The gds-proof engine already proves invariant preservation: given a predicate P(x) and a transition function x' = f(x, u), it checks whether P(x) ⟹ P(f(x, u)) via SymPy simplification. Lyapunov stability is exactly this pattern:

• Invariant 1: V(x) > 0 for all x ≠ 0 (positive definiteness)
• Invariant 2: V(f(x, u)) - V(x) < 0 (strict decrease along trajectories)

The proof machinery exists; we just need convenience constructors that translate control-theoretic specifications into the invariant format the engine already handles.

Proposed API

New file: gds_proof/analysis/lyapunov.py

from dataclasses import dataclass
import sympy

@dataclass(frozen=True)
class LyapunovResult:
    """Result of a Lyapunov stability analysis."""
    candidate: sympy.Expr              # V(x) expression
    positive_definite: str             # "PROVED" | "FAILED" | "INCONCLUSIVE"
    decreasing: str                    # "PROVED" | "FAILED" | "INCONCLUSIVE"
    stable: bool                       # True only if both proved
    dV_expr: sympy.Expr | None         # dV/dt or ΔV expression (for inspection)
    details: list[str]                 # human-readable proof trace

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
        where f(x) is given by state_transition.
    
    For discrete-time (continuous=False):
        Checks V(x) > 0 and V(f(x)) - V(x) < 0
        where f(x) is given by state_transition.
    
    Uses the existing gds-proof SymPy simplification strategies
    to attempt the proof.
    """

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 (eigenvalue check via SymPy)
    2. AᵀP + PA is negative definite (Lyapunov equation)
    
    This is the standard linear stability certificate.
    If P is the identity, this reduces to checking that A + Aᵀ
    is negative definite (sufficient but not necessary for stability).
    """

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).
    Solves symbolically via SymPy's linear system solver.
    Returns (P, result) if found, None if system is not stable.
    """

@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: str            # "PROVED" | "FAILED" | "INCONCLUSIVE"
    passive: bool
    details: list[str]

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.
    """

Implementation Notes

Integration with existing proof engine

The Lyapunov templates do NOT create new proof machinery. They translate to the existing patterns:

# lyapunov_candidate() internally does:
V_sym = parse_expr(V_expr)
# Construct the decrease condition as an invariant
if continuous:
    dVdt = sum(sympy.diff(V_sym, xi) * fi for xi, fi in zip(state_symbols, f_exprs))
    decrease_inv = dVdt < 0
else:
    V_next = V_sym.subs({xi: fi for xi, fi in zip(state_symbols, f_exprs)})
    decrease_inv = V_next - V_sym < 0

# Delegate to existing proof strategies
from gds_proof.analysis.symbolic import attempt_proof
pos_def_result = attempt_proof(V_sym > 0, assumptions=...)
decrease_result = attempt_proof(decrease_inv, assumptions=...)

The five existing strategies (VACUITY, DIRECT_SIMPLIFICATION, PREDICATE_IMPLICATION, Q_SYSTEM, BARE_IMPLICATION) handle the actual SymPy verification.

Limitations of symbolic approach

SymPy can prove stability for:

• Linear systems with small state dimension (2-4 states)
• Polynomial Lyapunov functions with simple dynamics
• Quadratic forms xᵀPx where P is given explicitly

SymPy will likely return INCONCLUSIVE for:

• High-dimensional systems (symbolic eigenvalue computation is O(n!) worst case)
• Non-polynomial dynamics (trig, exponential)
• Systems requiring SOS (sum-of-squares) decomposition

This is acceptable — the proof engine already returns INCONCLUSIVE as a valid result. Users can fall back to numerical eigenvalue checks in gds-analysis/linear.py for those cases.

Quadratic Lyapunov equation solver

find_quadratic_lyapunov() solves AᵀP + PA = -Q by:

1. Parameterizing P as a symmetric matrix with n(n+1)/2 unknowns
2. Expanding AᵀP + PA + Q = 0 into n² linear equations
3. Solving the linear system via sympy.solve()
4. Verifying the solution P is positive definite (eigenvalues > 0)

For small systems (n ≤ 4) this is fast. For larger systems, recommend using scipy's solve_continuous_lyapunov() in gds-analysis/linear.py instead.

Key Files

• packages/gds-proof/gds_proof/analysis/symbolic.py — five proof strategies (to be reused)
• packages/gds-proof/gds_proof/analysis/inductive_safety.py — three-layer inductive analysis (structural reference)
• packages/gds-proof/gds_proof/protocols.py — SymbolicBlock, SymbolicModel protocols
• packages/gds-proof/gds_proof/adapter.py — GDSSymbolicBlock adapter (reference for bridging)
• packages/gds-domains/gds_domains/symbolic/linearize.py — LinearizedSystem (A matrix source)

Testing

• test_lyapunov.py:
  • Stable linear: A = [[-1, 0], [0, -2]], V = x₁² + x₂² → PROVED stable
  • Unstable linear: A = [[1, 0], [0, -1]], same V → decrease FAILED
  • Quadratic: find_quadratic_lyapunov(A) for known stable A → returns valid P
  • Quadratic: unstable A → returns None
  • Nonlinear (Van der Pol near origin): V = x² + y², f = [y, -x + μ(1-x²)y] at μ=0 → PROVED (linear center)
  • Passivity: simple spring-mass with damping, V = ½(kx² + mv²), supply rate u·v → dissipation PROVED
  • INCONCLUSIVE case: high-dimensional system that exceeds SymPy's simplification capability → verify graceful result

Concepts Addressed (MATLAB Tech Talks)

• Video 1: Stability analysis (Lyapunov as alternative to eigenvalue methods)
• Video 8: Passivity-based control — energy dissipation, storage functions, supply rates, passivity theorem