Wide-Array of Nonlinear Dynamics Approximation

A zero-optimization, online streaming adaptive observer for real-time identification of nonlinear dynamical systems with time-varying parameters.

---

Overview

pywynda is a Python implementation of the WyNDA (Wide-Array of Nonlinear Dynamics Approximation) framework, a simultaneous state-and-parameter observer derived from the theory published in:

Hasan, A. (2024). WyNDA: A Method to Discover Mathematical Models of Nonlinear Dynamical Systems. MethodsX, 12, 102625. https://doi.org/10.1016/j.mex.2024.102625

The framework solves the nonlinear system identification problem in a strictly online, recursive fashion — consuming one measurement sample at a time with no batch history, no gradient descent, and no optimization solver of any kind. Its core engine is a dual-gain Kalman-type recursion that simultaneously corrects state estimates and parameter estimates at every time step, with exponential forgetting factors that make the observer adaptive to structural changes in the system dynamics.

---

Three Core Design Constraints

1. Zero-Optimization Constraint

No scikit-learn. No scipy.optimize. No neural network training loop. No gradient.

Every matrix operation reduces to a single numpy.linalg.solve call per time step. The observer is a pure matrix-algebra recursion — not a regression, not a minimization, and not a learned model. This makes it:

• Deterministic and reproducible without random seeds
• O(n³) per step in parameter dimension, with no iterative convergence
• Deployable on hardware that cannot run an optimization runtime

2. Strict Online Streaming

One sample in → corrected state + updated parameters out.

The .step(y, Psi) method maintains all internal state between calls. There is no buffer, no window, no queue. The observer is fully compatible with real-time sensor streams, embedded control loops, and hardware-in-the-loop pipelines.

3. Non-Stationary Adaptation Edge

Exponential forgetting discounts stale history. No retraining required.

Two scalar forgetting factors — $\lambda_x \in (0,1)$ for state estimation and $\lambda_\theta \in (0,1)$ for parameter estimation — inflate the respective covariance matrices at each step, ensuring that recent data always dominates the gain computation. When system parameters change abruptly (component failure, structural drift, regime transition), the observer automatically re-identifies the new dynamics within a time window of approximately $\tau \approx \Delta t ,/,(1 - \lambda_\theta)$.

---

Mathematical Engine

The observer solves the following discrete-time system identification problem at each step $k$.

System Model

Equation	Role
$y(k) = x(k) + v(k)$	Measurement with noise $v(k)$
$x(k+1) = x(k) + \Psi(y(k), u(k)),\theta$	Euler-discretised nonlinear dynamics

where $\Psi \in \mathbb{R}^{n \times p}$ is the block-diagonal regressor library evaluated at the current measurement $y(k)$ and control input $u(k)$, and $\theta \in \mathbb{R}^p$ is the unknown parameter vector.

Gain Computations

$$K_{x}(k) = P_{x}(k|k-1)\bigl(P_{x}(k|k-1) + R_{x}(k)\bigr)^{-1}$$

$$\Omega(k) = \Gamma(k|k-1),P_{\theta}(k|k-1),\Gamma^{\top}(k|k-1) + R_{\theta}(k)$$

$$K_{\theta}(k) = P_{\theta}(k|k-1),\Gamma^{\top}(k|k-1),\Omega(k)^{-1}$$

where $\Gamma(k|k-1)$ is the coupling matrix that links the parameter subspace to the innovation signal, and $R_x$, $R_\theta$ are measurement noise covariances acting as natural Tikhonov-style regularisers.

State and Parameter Correction

$$\bar{x}(k|k) = \bar{x}(k|k-1) + \bigl(K_{x}(k) + \Gamma(k|k-1),K_{\theta}(k)\bigr)\bigl(y(k) - \bar{x}(k|k-1)\bigr)$$

$$\bar{\theta}(k|k) = \bar{\theta}(k|k-1) - K_{\theta}(k)\bigl(y(k) - \bar{x}(k|k-1)\bigr)$$

One-Step-Ahead Predictions and Covariance Transitions

$$\bar{x}(k+1|k) = \bar{x}(k|k) + \Psi(y(k), u(k)),\bar{\theta}(k|k)$$

$$P_{x}(k+1|k) = \lambda_{x}^{-1}\bigl(I - K_{x}(k)\bigr)P_{x}(k|k-1)$$

$$P_{\theta}(k+1|k) = \lambda_{\theta}^{-1}\bigl(I - K_{\theta}(k),\Gamma(k|k-1)\bigr)P_{\theta}(k|k-1)$$

$$\Gamma(k+1|k) = \bigl(I - K_{x}(k)\bigr)\Gamma(k|k-1) - \Psi(y(k), u(k))$$

The $\lambda^{-1}$ inflation in the covariance transitions is the mechanism of exponential forgetting — each step slightly enlarges $P_x$ and $P_\theta$, increasing the observer gain and responsiveness to new data.

Theorem 1 (MethodsX 2024): Exponential Error Decay

Under the Persistency of Excitation (PE) condition — that the regressor $\Psi$ spans the full parameter space over a finite time window — both error norms decay exponentially:

$$|\varepsilon_x(k)| \leq C_x,e^{-\alpha_x k}, \qquad |\varepsilon_\theta(k)| \leq C_\theta,e^{-\alpha_\theta k}$$

This guarantee is validated empirically across all five benchmark examples in this repository.

---

Repository Layout

pywynda/
├── src/
│   └── pywynda/
│       ├── __init__.py          # Public API surface
│       ├── core.py              # WyNDAObserver — the recursive gain engine
│       ├── utils.py             # Trajectory generation, PE checks, error metrics
│       └── py.typed             # PEP 561 marker (fully typed)
│
├── tests/
│   └── test_observer.py         # 26 unit tests across 8 validation categories
│
├── examples/
│   ├── ex1_msd.py               # Linear MSD — over-parameterised library
│   ├── ex2_nonlinear.py         # Lorenz, Rossler, Lotka-Volterra, Van der Pol
│   ├── ex3_control_pe.py        # PE vs. no-PE comparison on forced oscillator
│   ├── ex4_maglev.py            # Magnetic levitation — 3-state physical system
│   ├── ex5_sindy_comparison.py  # WyNDA vs. Batch SINDy: non-stationary drift
│   ├── run_all.py               # Sequential manifest runner for all examples
│   └── plots/                   # Auto-generated PNG convergence/comparison plots
│
├── pyproject.toml
├── uv.lock
└── README.md

---

Installation

pywynda uses the uv toolchain for environment management.

# Clone the repository
git clone https://github.com/scarwizz/pywynda.git
cd pywynda

# Create and activate the virtual environment
uv venv
# Windows:
.venv\Scripts\activate
# Linux / macOS:
source .venv/bin/activate

# Install all dependencies (numpy, scipy, matplotlib, pytest)
uv sync

# Verify installation with the full test suite
uv run pytest tests/ -v

Runtime dependencies: numpy >= 2.5, scipy >= 1.18, matplotlib >= 3.8 No scikit-learn, torch, tensorflow, jax, or any optimization framework.

---

Quickstart

import numpy as np
from pywynda import WyNDAObserver

# ── System dimensions ──────────────────────────────────────────────────────
n_states = 2   # [position, velocity]
n_params = 6   # 3 basis functions per state equation (block-diagonal library)
dt       = 0.001

# ── Initialise the observer ────────────────────────────────────────────────
obs = WyNDAObserver(
    n_states   = n_states,
    n_params   = n_params,
    lambda_x   = 0.999,          # state-estimation forgetting factor
    lambda_theta = 0.998,        # parameter-estimation forgetting factor (tau ~ 0.5 s)
    P_x0       = 100.0 * np.eye(n_states),   # initial state covariance
    P_theta0   = 500.0 * np.eye(n_params),   # initial parameter covariance
    R_x        = 0.05**2 * np.eye(n_states), # measurement noise covariance
    R_theta    = 5e-3  * np.eye(n_states),   # parameter noise covariance (n x n)
    x0_hat     = np.zeros(n_states),         # initial state estimate
    theta0_hat = np.zeros(n_params),         # initial parameter guess
)

# ── Build the block-diagonal regressor Psi (n x p) at measurement y ───────
def build_psi(y: np.ndarray) -> np.ndarray:
    x1, x2 = y[0], y[1]
    phi = dt * np.array([x1, x2, x1**2])     # 3-entry basis per equation
    Psi = np.zeros((n_states, n_params))
    Psi[0, :3] = phi                          # x1 equation
    Psi[1, 3:] = phi                          # x2 equation
    return Psi

# ── Streaming loop — one sample at a time ─────────────────────────────────
# (replace with your real sensor stream)
for k, y_k in enumerate(your_measurement_sequence):
    state = obs.step(y=y_k, Psi=build_psi(y_k))

    x_hat     = state.x_hat      # corrected state estimate, shape (n,)
    theta_hat = state.theta_hat  # updated parameter vector, shape (p,)
    inno      = state.innovation # innovation signal y(k) - x_hat(k|k-1)

# ── Known-physics feed-forward (optional) ─────────────────────────────────
# Inject analytically-known contributions (e.g. gravity, kinematic coupling)
# without contaminating the parametric subspace:
obs.inject_feedforward(delta_x=np.array([dt * float(y_k[1]), dt * 9.81]))

---

Benchmark Examples

Run all five examples sequentially and regenerate all plots:

uv run python examples/run_all.py

Or run individual examples:

uv run python examples/ex1_msd.py          # ~3 s
uv run python examples/ex2_nonlinear.py    # ~25 s (4 chaotic systems)
uv run python examples/ex3_control_pe.py   # ~5 s
uv run python examples/ex4_maglev.py       # ~3 s
uv run python examples/ex5_sindy_comparison.py  # ~10 s

Example Suite Summary

Example	System	Key Result
Ex 1	Mass-Spring-Damper (over-parameterised, 12 basis)	$k/m$ and $b/m$ recovered; spurious coefficients collapse toward zero
Ex 2a	Lorenz strange attractor	$\hat{\sigma}=9.78$, $\hat{\rho}=27.79$, $\hat{\beta}=2.98$; $\varepsilon_\theta = 0.53$ at $t=2$ s
Ex 2b	Rössler system	Expected PE-deficit failure on $c$-parameter — documented in §4.2 of reference
Ex 2c	Lotka-Volterra predator-prey	All 4 parameters within 5 % of truth; $\varepsilon_\theta(\text{final}) = 0.034$
Ex 2d	Van der Pol oscillator	$\hat{\mu} = 1.196$ (0.3 % error); dual-coefficient consistency confirmed
Ex 3	Forced oscillator — PE vs. no-PE	Empirical proof of Theorem 1 necessity condition
Ex 4	Magnetic levitation (3-state)	$\hat{L}=3.002$ H, $\hat{R}=6.004\ \Omega$, $\hat{C}=1.497$ F (sub-1 % error in 0.2 s)
Ex 5	MSD with sudden parameter drift	WyNDA vs. Batch SINDy — see dedicated section below

---

Benchmark Showcase: WyNDA vs. Batch SINDy under Non-Stationary Dynamics

Example 5 (ex5_sindy_comparison.py) is a rigorous head-to-head evaluation on the most challenging scenario in data-driven system identification: sudden parameter drift mid-operation.

Experiment Design

A Mass-Spring-Damper system is subjected to a hard step change in spring stiffness at $t = 2.5$ s, simulating a sudden structural failure or component degradation:

$$k(t) = \begin{cases} 84.0 \text{ N/m} & 0.0 \leq t < 2.5 \text{ s} \ 40.0 \text{ N/m} & 2.5 \leq t \leq 5.0 \text{ s} \end{cases}$$

Multi-frequency PE forcing ($3\sin(2\pi \cdot 1.5, t) + 2\sin(2\pi \cdot 3.7, t)$) ensures the PE lower bound is satisfied throughout. White Gaussian noise ($\sigma = 0.05$) is applied to all state measurements.

Both algorithms run on the exact same trajectory (identical random seed):

Parameter	WyNDA	Batch SINDy
Processing mode	Online streaming (1 sample at a time)	Batch (re-fit every 200 steps over all history)
Forgetting	$\lambda_\theta = 0.998$ ($\tau \approx 0.5$ s)	None — accumulates all data from $t=0$
Solver	np.linalg.solve (one call/step)	np.linalg.solve on $\Theta^\top\Theta + \alpha I$
Sklearn / optimize	No	No (zero-optimization compliant)

Performance Telemetry

Metric	WyNDA	Batch SINDy	Advantage
Full-run $k$ RMSE	20.4 N/m	35.3 N/m	WyNDA 1.7× lower
Post-drift $k$ RMSE	23.3 N/m	41.3 N/m	WyNDA 1.8× lower
$\hat{k}$ at $t = 5.0$ s	41.5 N/m ✓	79.4 N/m ✗	WyNDA converges; SINDy does not
Detection latency (5 N/m tol.)	1040 ms	> T_end (never)	SINDy structurally fails

Why SINDy Fails Here

SINDy's batch regression accumulates all historical data $[0 \ldots k]$ into the design matrix $\Theta$. After the drift at $t = 2.5$ s, the regression is contaminated by 2 500 steps of pre-drift dynamics where $k = 84$ N/m. The normal equations mix two incompatible physical regimes and converge to a weighted mean of the two spring constants — stuck near 79 N/m at end-of-run — regardless of how many post-drift samples arrive.

WyNDA's forgetting factor $\lambda_\theta = 0.998$ means data from $3\tau = 1.5$ s ago carries only $\lambda_\theta^{1500} \approx 5%$ weight in the gain computation. The observer effectively erases the pre-drift regime within two forgetting time-constants and re-identifies $k = 40$ N/m cleanly.

This demonstrates WyNDA's definitive architectural edge: global batch regression is structurally incompatible with non-stationary dynamical systems. Online streaming with forgetting is not a heuristic optimisation — it is the mathematically correct formulation for this class of problems.

---

Test Suite

The full validation suite confirms all theoretical guarantees from MethodsX 12 (2024) 102625.

========================= test session starts ==========================
platform win32 -- Python 3.13.11, pytest-9.1.1, pluggy-1.6.0
rootdir: C:\...\pywynda
configfile: pyproject.toml
collected 26 items

tests/test_observer.py::TestOutputShapes::test_observer_state_shapes         PASSED [  3%]
tests/test_observer.py::TestOutputShapes::test_observer_history_shapes        PASSED [  7%]
tests/test_observer.py::TestOutputShapes::test_psi_sequence_shape             PASSED [ 11%]
tests/test_observer.py::TestConvergence::test_state_estimation_convergence    PASSED [ 15%]
tests/test_observer.py::TestConvergence::test_parameter_estimation_convergence PASSED [ 19%]
tests/test_observer.py::TestConvergence::test_final_parameter_accuracy         PASSED [ 23%]
tests/test_observer.py::TestConvergence::test_state_error_last_quarter_small  PASSED [ 26%]
tests/test_observer.py::TestConvergence::test_exponential_decay_envelope      PASSED [ 30%]
tests/test_observer.py::TestPersistenceOfExcitation::test_chirp_satisfies_pe  PASSED [ 34%]
tests/test_observer.py::TestPersistenceOfExcitation::test_white_noise_satisfies_pe PASSED [ 38%]
tests/test_observer.py::TestPersistenceOfExcitation::test_gram_matrix_positive_definite PASSED [ 42%]
tests/test_observer.py::TestDeterminism::test_same_seed_same_results          PASSED [ 46%]
tests/test_observer.py::TestDeterminism::test_different_seeds_different_results PASSED [ 50%]
tests/test_observer.py::TestForgettingFactor::test_lower_lambda_faster_initial_theta_convergence PASSED [ 53%]
tests/test_observer.py::TestPredictionIntegrity::test_x_pred_finite           PASSED [ 57%]
tests/test_observer.py::TestPredictionIntegrity::test_theta_pred_finite       PASSED [ 61%]
tests/test_observer.py::TestPredictionIntegrity::test_innovations_bounded     PASSED [ 65%]
tests/test_observer.py::TestPredictionIntegrity::test_covariances_decrease    PASSED [ 69%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_x_zero       PASSED [ 73%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_x_too_large  PASSED [ 76%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_theta_negative PASSED [ 80%]
tests/test_observer.py::TestInputValidation::test_wrong_P_x_shape             PASSED [ 84%]
tests/test_observer.py::TestInputValidation::test_run_mismatched_steps        PASSED [ 88%]
tests/test_observer.py::TestRunAPI::test_run_produces_all_fields              PASSED [ 92%]
tests/test_observer.py::TestRunAPI::test_run_step_consistency                 PASSED [ 96%]
tests/test_observer.py::TestRunAPI::test_multiple_input_types_converge        PASSED [100%]

========================= 26 passed in 4.77s ===========================

Test Category Coverage

#	Category	Scenarios	Validates
1	Shape Tests	3	All output matrices carry correct dimensions
2	Convergence Tests	5	Exponential decay of $\varepsilon_x$ and $\varepsilon_\theta$ (Theorem 1)
3	Persistency of Excitation	3	PE lower bound satisfied for chirp and white-noise inputs
4	Determinism	2	Same seed → bit-identical trajectories
5	Forgetting Factor Sensitivity	1	Lower $\lambda$ gives faster initial convergence
6	Prediction Integrity	4	$\hat{x}(k+1
7	Input Validation	5	Constructor and .step() raise on bad arguments
8	Run() API	3	ObserverHistory consistency with step-by-step calls

---

Public API Reference

WyNDAObserver

class WyNDAObserver:
    def __init__(
        self,
        n_states:     int,
        n_params:     int,
        lambda_x:     float,          # state forgetting factor, in (0, 1)
        lambda_theta: float,          # parameter forgetting factor, in (0, 1)
        P_x0:         NDArray,        # initial state covariance, shape (n, n)
        P_theta0:     NDArray,        # initial parameter covariance, shape (p, p)
        R_x:          NDArray,        # measurement noise covariance, shape (n, n)
        R_theta:      NDArray,        # parameter noise covariance, shape (n, n)
        x0_hat:       NDArray,        # initial state estimate, shape (n,)
        theta0_hat:   NDArray,        # initial parameter estimate, shape (p,)
    ) -> None: ...

    def step(self, y: NDArray, Psi: NDArray) -> ObserverState:
        """Process one measurement. Returns the corrected state and parameters."""

    def inject_feedforward(self, delta_x: NDArray) -> None:
        """Add a known non-parametric dynamics correction to the one-step-ahead prior.
        Use for gravity, kinematic coupling, or other analytically-known contributions."""

    def run(self, Y: NDArray, Psi_seq: NDArray) -> ObserverHistory:
        """Process a full measurement sequence and return stacked history."""

ObserverState (returned by .step())

Field	Shape	Description
x_hat	(n,)	Corrected state estimate $\bar{x}(k
theta_hat	(p,)	Updated parameter vector $\bar{\theta}(k
x_pred	(n,)	One-step-ahead prediction $\bar{x}(k+1
theta_pred	(p,)	One-step-ahead parameter prediction
P_x	(n,n)	State covariance $P_x(k+1
P_theta	(p,p)	Parameter covariance $P_\theta(k+1
Gamma	(n,p)	Coupling matrix $\Gamma(k+1
K_x	(n,n)	State Kalman gain $K_x(k)$
K_theta	(p,n)	Parameter Kalman gain $K_\theta(k)$
innovation	(n,)	Innovation signal $y(k) - \bar{x}(k

---

Design Principles

Regressor Library Psi — Block-Diagonal Structure

The regressor $\Psi \in \mathbb{R}^{n \times p}$ must be block-diagonal to decouple each state equation's parameter subspace. For an $n$-state system with $m$ basis functions per equation:

Psi = [ phi(y, u) |   0     ]   Row 0: x1-equation parameters
      [   0       | phi(y,u)]   Row 1: x2-equation parameters
                   p = n * m

This structure ensures the parameter correction $K_\theta \cdot \text{innovation}$ updates only the relevant equation coefficients at each step.

Forgetting Factor Selection Guide

Regime	$\lambda_\theta$ recommendation	Effective memory window
Stationary system	$0.999$	~1 000 steps
Slow drift	$0.998$	~500 steps
Sudden step-change	$0.995$ – $0.998$	200–500 steps
Rapid adaptation needed	$0.990$	~100 steps

Setting $\lambda_\theta$ too low causes covariance wind-up during stable phases; too high causes slow drift tracking. The sweet spot for step-change detection is $\lambda_\theta \approx 0.998$, giving a forgetting time-constant $\tau = \Delta t / (1 - \lambda_\theta)$.

Solver Compliance

Every linear system in the recursion is solved via numpy.linalg.solve, never numpy.linalg.inv. This follows numerical best practice: solving $A x = b$ directly via LU factorisation is both faster and numerically stabler than computing $A^{-1} b$.

---

Reference

@article{hasan2024wynda,
  title   = {{WyNDA}: A Method to Discover Mathematical Models of
              Nonlinear Dynamical Systems},
  author  = {Hasan, Agus},
  journal = {MethodsX},
  volume  = {12},
  pages   = {102625},
  year    = {2024},
  doi     = {10.1016/j.mex.2024.102625},
  url     = {https://doi.org/10.1016/j.mex.2024.102625}
}

---

License

MIT License — see LICENSE for details.

---

_{Built with NumPy. Zero optimization dependencies. Strictly online. Reference implementation of MethodsX 12 (2024) 102625.}