# Online Stability Indicator (OSI) — Proposed Solutions

> Folder note: the requested folder name is "claude trial"; it is created as
> `claude_trial` (underscore) so the path is usable as a MATLAB folder and
> consistent with the rest of the repository.

## 1. Problem restated

We want an **online stability indicator** computed from streaming measurements
that detects shrinkage of the closed-loop **NCF (normalized coprime factor)
stability margin** during an experiment, so an operator can react (switch
controller, pause) *before* the loop goes unstable.

### Loop topology (from `sys_id.m`)

```
 dd ──►(+)──► e ──► C(z) ──► dc ──► H_dc_dm ──► dm ──┐
        ▲ -                  (PID)     (plant)        │
        └────────────────────────────────────────────┘   (unity feedback)
```

* `dd` — desired displacement (reference, exogenous excitation)
* `e  = dd − dm` — loop error (open-loop input)
* `dc` — digital-controller output (command)
* `dm` — measured displacement
* `C(z)` — discrete PID, **constant and known**: `P=0.15, I=20, D=0.0015,
  N=512`, Forward-Euler, `Ts = 1/2048 s`
* `H_dc_dm` — actuator + specimen plant; **drifts over the experiment**

Open loop: `L = C·H_dc_dm`. Closed-loop relations used by the estimators:
`dm = T·dd`, `e = S·dd`, with `S = 1/(1+L)`, `T = L/(1+L)`, and `dm = L·e`.

### Margin definition (reuse of `custom_stability_margins.m`)

The repository's loop-only NCF margin is the minimum normalized Nyquist
distance of the loop over a frequency band:

```
b_looponly = min_ω  |1 + L(jω)| / sqrt( 2·(1 + |L(jω)|²) )
```

This is the value all estimators below reproduce online. (The rigorous
4‑block Vinnicombe/Glover–McFarlane margin `b_{P,C}` from Robust Control
Toolbox `ncfmargin` is also reported by Methods 2 and 4.)

### Facts established from the data (see `report.md` for the numbers)

* `fs = 2048 Hz`; test record 823.5 s; analysis window **620 → 770.1 s**.
* Nominal loop-only NCF = **0.588**, and it is **set at ≈ 35 Hz** (a lightly
  damped actuator/specimen resonance), *not* at the 1.36 Hz gain crossover.
  → the window/FFT must resolve ≈ 35 Hz, which it easily does.
* `dd` has broadband energy (50 % below 41 Hz), so ≈ 35 Hz is well excited.
* At the start of the window the measured loop matches the `exp_tf` model to
  within ≈ 5 % → the "model correct at the beginning" assumption holds.
* `rms(dm)` is steady (~0.08) until ≈ 700 s, then climbs and triples by
  760–770 s → the margin shrinkage we must detect (true instability/divergence
  is just past 770 s, where `dm` reaches −14.6).

### Windowing

* Each indicator sample uses one **non-overlapping** window of length `W` (s)
  tiling 620→770.1 s (`W = 10 s` default → 15 samples). Non-overlap makes the
  indicator samples (approximately) independent — the white-measurement
  assumption a downstream **Kalman filter** needs. Internal Welch averaging may
  still overlap *inside* a window.

---

## 2. The five proposed methods

All methods output the same deliverable: a **time series of the estimated NCF
margin** over 620→770.1 s. They differ in how the loop `L(jω)` (or the plant
`H_dc_dm`) is estimated inside each window.

### Method 1 — Nonparametric spectral, direct (H1)  →  `trial_1`
Per window, estimate the plant FRF by the H1 estimator from the *measured*
command and response,
`Ĥ(jω) = P_{dc,dm}(ω) / P_{dc,dc}(ω)` (Welch, Hann),
form `L̂ = C(jω)·Ĥ(jω)` with the **known** PID, and take the loop-only NCF.
* **Pros:** model-free, fast, no tuning, faithful to the resonance (gain
  weighting suppresses the noisy high-frequency tail).
* **Cons:** `dc` is generated by feedback, so H1 is mildly **closed-loop
  biased** (output noise leaks into the input) → small variance in the stable
  region.
* **Role:** simple, dependable baseline.

### Method 2 — Parametric per-window model ID  →  `trial_2`
Per window, fit a continuous rational plant `Ĥ(s)` (order 5p/3z, matching the
`exp_tf` structure) to the Welch FRF with `tfest` (`EnforceStability`), exactly
as `sys_id.m` does, then `c2d` to `Ĥ(z)` and `L̂ = C·Ĥ`. Reports **both** the
loop-only NCF and the rigorous `ncfmargin` `b_{P,C}`, plus closed-loop poles.
* **Pros:** smooth rational loop; gives poles and the exact 4-block margin;
  physically interpretable.
* **Cons:** per-window batch ID on 10 s of noisy closed-loop data is
  **higher-variance / occasionally over-optimistic** (a poor fit smooths the
  resonance dip and *over*states the margin — dangerous for a safety monitor);
  slowest; occasional bad fits.

### Method 3 — Closed-loop indirect, instrument-variable (dd)  →  `trial_3`
Use the exogenous reference `dd` as an instrument to remove closed-loop bias:
`Ĥ(jω) = P_{dd,dm}(ω) / P_{dd,dc}(ω)` (ratio of cross-spectra), then
`L̂ = C·Ĥ`, loop-only NCF.
* **Pros:** statistically the **right** estimator for closed-loop data →
  lowest variance, starts essentially on the nominal value; still model-free.
* **Cons:** relies on `dd` being persistently exciting and uncorrelated with
  plant noise; degrades if the reference content collapses (e.g., once the
  loop self-oscillates).

### Method 4 — Model-reference ν-gap drift / early-warning  →  `trial_4`
The only method that uses the **`exp_tf` model as the baseline** (the
objective's "model correct at the beginning" assumption). Per window fit a
5p/3z plant, then measure the **Vinnicombe ν-gap** to the window-1 baseline
(self-referenced to cancel common estimation bias), `δ_k = δ_ν(Ĥ_1, Ĥ_k)`, and
report the robust-stability prediction `b̂_k = b_{P,Ĥ_1} − δ_k` (median-of-3
robustified). Alarm when `δ_k → b_{P,Ĥ_1}`.
* **Pros:** directly encodes the robust-stability theorem (loop stays stable
  while `δ_ν < b_{P,C}`); gives the **earliest** warning; a clean
  switch/pause trigger.
* **Cons:** ν-gap integrates model mismatch over all frequencies → **most
  fragile** on short noisy windows; better as a binary early-warning than a
  precise margin; needs smoothing.

### Method 5 — Recursive online tracking (RLS / Kalman)  →  `trial_5`
Run a recursive-least-squares ARX tracker (`na=nb=5`, forgetting factor `λ`) on
decimated (256 Hz) `dc→dm` data continuously; at each non-overlapping window
boundary snapshot the parameters, build `Ĥ(z)`, `L̂ = C·Ĥ`, loop-only NCF.
* **Pros:** truly **online / O(1) per sample**, naturally smooth, lowest
  latency between samples, and the most "Kalman-ready" (the parameter tracker
  *is* a KF on a random-walk model). Decimation conditions the estimation.
* **Cons:** forgetting factor introduces **lag** and the recursion inherits
  closed-loop bias → tends to read **optimistically high** and detect a little
  late; `λ` must be tuned (responsiveness vs. smoothness).

### Method 6 — Recursive frequency-point Kalman filter (fast, small window)  →  `trial_6`
*(added to answer "make the window much smaller / detect within 0.5 s".)*
Place one complex Kalman state on each of a few discrete TF points `H(f_i)` near
the binding ~35 Hz resonance, initialised from the **known model**. Every
**0.5 s**, form a short-window H1 measurement `z_i` and coherence `γ_i²` at each
point, and fuse it with a Kalman update whose **measurement noise R_i is adaptive
in the coherence** (`R_i ∝ (1−γ_i²)/(n_d γ_i²)·|z_i|²`): well-excited points are
trusted, poorly-excited ones ignored. The margin is a robust percentile of the
loop distance `L_i = C(f_i)H_i`.
* **Pros:** **0.5 s time resolution** (20× M1–M5) and update latency; reliability
  comes from *temporal* Kalman fusion, not in-window averaging, so the window can
  be tiny; smoothest stable phase of all; uniquely **resolves rapid collapse**.
* **Cons:** short buffers under-resolve the sharp 35 Hz dip (leakage) → slightly
  optimistic on slow erosion; best as a **fast trip** beside M3's accurate slow
  estimate. Full derivation in `method_6_explained.md`.

### Method 7 — Output-only closed-loop resonance / sensitivity-peak margin  →  `trial_7`
*(added to catch a 35 Hz instability under a **low-frequency input** — the 1749
file — where the open-loop FRF at 35 Hz is a tautology and M1–M6 are blind.)*
The binding 35 Hz mode is excited *from inside the loop* (its growth **is** the
instability), so it is visible in the response `dm` **alone**. Estimate the
closed-loop resonance amplification `Spk = √(peak/floor PSD of dm)` ≈ peak
sensitivity `‖S‖∞ ≈ 1/min|1+L|`, then map to the loop-only NCF, **calibrated to
the model's nominal margin at the start**: `b(t) = b_nom · Spk_ref / Spk(t)`.
* **Pros:** **input-agnostic** — works for broadband (1731) and low-frequency
  (1749) inputs; catches the instability from the output only; cross-validates
  the FRF methods on 1731 and is the *only* method that catches 1749.
* **Cons:** a *relative* (start-calibrated) margin, not an absolute one; needs a
  longer (4 s) buffer for a stable resonance estimate; assumes the binding mode
  shows as a resonance in `dm`. Detailed in `common/osi_method7_resonance.m`.

### Cross-cutting — Kalman post-smoothing
Because the windows are non-overlapping, the indicator series is a clean
measurement sequence `z_k`. A scalar random-walk Kalman filter (state = true
margin, `m_k = m_{k-1} + w_k`, `z_k = m_k + v_k`) — optionally with an RTS
smoother — is applied to the recommended method's output to denoise the
indicator and provide a one-step-ahead prediction. This is demonstrated in
`run_all_and_compare.m`.

---

## 3. How each trial is validated as "complete"

For every method the trial script:
1. loads the test record and the `exp_tf` nominal model + PID,
2. computes the margin for every non-overlapping window in 620→770.1 s,
3. **saves** `results/trial_k_result.mat` and a PNG of the series,
4. prints the per-window table.

A method "achieves the expected outcome" iff it returns a full margin time
series that (a) starts near the nominal 0.588 (model correct at start) and
(b) declines as the loop degrades toward 770 s. `run_all_and_compare.m`
verifies this for all five and quantifies accuracy, smoothness, and detection
latency.