SOL‑OSI: a loop that measures its own distance to instability

01 Setting

The experiment that must not break

In a real-time hybrid simulation (RTHS), half of a structure exists as software and half as hardware. A finite-element model of a building computes, in real time, how a floor should move during an earthquake. A hydraulic actuator then imposes that motion on a physical column specimen in the lab, measures the force the column pushes back with, and feeds that force into the model. Software and steel take turns, thousands of times per second.

Inside the machinery sits an ordinary servo loop. A controller receives a displacement command and drives the actuator so that the measured displacement tracks it. Two recorded signals are all SOL‑OSI ever sees:

Here is the danger. The column specimen changes during the test — that is the point of testing it. Damage shifts its stiffness and damping, which shifts the loop’s dynamics, which erodes its stability margin. If the loop crosses into instability, the actuator oscillates violently and the specimen — possibly months of preparation — is destroyed by the test equipment rather than tested by it.

So we want an online stability indicator (OSI): a number b, recomputed every 0.125 s from u and y alone, that measures how far the loop currently is from instability. Large b: comfortable. Falling b: the loop is drifting toward the edge, and a supervisor should gently turn the experiment down before physics turns it off. SOL‑OSI is the estimator of b; the supervisor is the authority regulator we meet at the end.

1function [out, state, info] = sol_osi(varargin)
2%SOL_OSI Healthy-reference two-branch NCF-margin estimator.
3%   cfg = sol_osi('config', dt, initial_loop_tf, override)
4%   state = sol_osi('state', cfg)
5%   [b, state, info] = sol_osi(dd, dm, cfg, state)

02 The model

The loop on paper

Start from what you know. Each block in the diagram is a transfer function: the ratio of output to input in the Laplace domain, G(s) = Y(s)/U(s). Evaluating it on the imaginary axis, s = jω, gives the frequency response — a complex number per frequency whose magnitude is gain and whose angle is phase shift.

The controller in this rig is the classic PID with a filtered derivative, with the exact constants used on the hardware:

The plant H(s) — hydraulic actuator, servo valve, and the column itself, from displacement command to measured displacement — is not written down from first principles. It was identified from measurements on the rig and stored in the codebase as H_dc_dm (the 7k column model). Multiplying the two blocks that the signal passes through gives the open-loop transfer function, the single object the whole module revolves around:

SOL-OSI never works with a formula for L₀. At configuration time it samples the frequency response on a fixed grid: the band 2–100 Hz in steps of 2 Hz — 50 complex numbers. Two sub-bands get special treatment: 58–62 Hz is masked out (mains electrical interference lives there, and a stability estimator should not be fooled by the power grid), and 28–44 Hz is reserved as the residual fit band for Section 07. Below is the actual magnitude data from the rig:

62f = (0:cfg.segment_n/2)'*cfg.fs/cfg.segment_n;
63cfg.bin = find(f >= cfg.band_hz(1) & f <= cfg.band_hz(2));
64cfg.freq_hz = f(cfg.bin);
65cfg.keep = ~(cfg.freq_hz >= cfg.line_reject_hz(1) & ...
66    cfg.freq_hz <= cfg.line_reject_hz(2));
72cfg.L0 = squeeze(freqresp(initial_loop_tf, 2*pi*cfg.freq_hz));

03 Two lines of algebra

Closing the loop

The loop is closed, so the signals we record are closed-loop signals. Before we can measure anything we need the map from command u to measurement y. Follow the block diagram around once — each line below is one substitution:

1. The error is what enters the controller: command minus measurement. e = u − y
2. The measurement is the error pushed through both blocks, C then H — that product is L. y = C H e = L e = L (u − y)
3. Collect the y terms on the left. (1 + L) y = L u
4. Divide. The closed-loop transfer function from command to measurement: T = yu = L1 + L
5. And the error itself obeys a twin formula — divide line 3’s logic into line 1: S = eu = 11 + L ,  S + T = 1

T is the complementary sensitivity — how faithfully y tracks u. S is the sensitivity — how much of the command leaks into the error, and, just as important for us, how much the loop amplifies disturbances injected at the output. Both formulas are evaluated frequency by frequency: plug in the complex number L(jω), get complex T(jω) and S(jω).

04 The central idea

Distance to instability

Section 03 ended with the observation that instability means L(jω) touching −1. The most natural margin, then, is simply a distance in the complex plane, computed at every grid frequency:

This is a fine start — it is exactly the distance from the Nyquist curve to the critical point, and its reciprocal is the peak of |S|. But it has a flaw for a monitor: it treats all parts of the plane unequally. Out where |L| is huge, the curve can swing by whole units of distance while the loop’s behavior barely changes; near the origin, a small step is a big deal. A margin that will be compared against one fixed threshold for the whole test needs a distance that weighs every neighborhood fairly.

The fix is to measure distance on a normalized version of the plane. Divide out the size of the points involved:

1. The chordal distance between two complex numbers z₁ and z₂ shrinks each leg by its own magnitude scale: c(z₁, z₂) = z₁ − z₂ 1 + z₁²1 + z₂²
2. Set z₁ = L(jω) and z₂ = −1. The second root becomes √2, a constant: c(L, −1) = 1 + L(jω) 2(1 + L(jω)²)
3. This is SOL-OSI’s per-frequency margin, capped at its own largest possible value: b(ω) = min( 1 + L 2 (1 + L²),  12 )

Two consequences of the formula are worth pausing on. First, the cap: when the loop is switched off entirely, L = 0, the distance to −1 is 1/2 ≈ 0.707 — no loop can ever be “safer” than no loop, so b lives in [0, 0.707]. Second, the loop is only as safe as its worst frequency, so the published margin takes the minimum across the grid, skipping the masked mains bins:

Everything after this section is bookkeeping in service of one question: what is the actual, current L(jω) to feed into this formula? Before that, get a feel for the geometry. The explorer below applies two knobs to the real loop: a gain a (the “authority” — Section 09 will turn this knob automatically) and a factor ρ standing for damage-induced loop-gain change:

293function b = chordal(L, cap)
294b = min(real(abs(1+L)./sqrt(2*(1+abs(L).^2))), cap);
295end
297function y = minMasked(x, mask)
298y = min(x(mask));
299end

05 Measurement

Listening to the loop

The formula for b wants L(jω) — the loop as it is now, damage included. That must come from the two live signals. SOL-OSI processes them in small blocks: every 0.125 s it receives 256 fresh samples and updates a sliding 0.5 s buffer (1024 samples) of each signal. The buffer, not the block, is what gets analyzed — a 0.5 s window is what sets the 2 Hz frequency grid you have already seen (Δf = 1/0.5 s = 2 Hz).

Each analysis starts by tapering the buffer with a Hann window and taking the FFT of both signals:

1. The window fades the buffer in and out so the FFT does not see false jumps at the buffer edges (spectral leakage): w[n] = 0.5 − 0.5 cos (2πnN − 1) ,  N = 1024
2. Transform the windowed signals; keep only the 2–100 Hz bins: U_k = FFT{w · u} |_k ,  Y_k = FFT{w · y} |_k
3. From the pair, form three real-time statistics per bin — two powers and one cross term: p_u = U_k² ,  p_y = Y_k² ,  c_yu = Y_k U_k^*

One 0.5 s window is far too noisy to trust. The classical remedy is to average many windows; the real-time remedy is an exponentially weighted moving average (EWMA) — a one-line filter that keeps a running mean which gently forgets the past:

Applied to the three statistics with τ = 0.25 s (so α = 1 − e^−0.5 = 0.393), this yields the smoothed auto-spectra and cross-spectrum that every later formula consumes:

One more scalar is computed per block, and it will gate almost everything downstream. The estimator asks: is there any excitation to learn from at all?

The fourth power makes this a soft switch: drive comfortably above the quiet threshold 0.06 and activity saturates at 1; drive below it and activity plunges toward 0. When the room is quiet, spectral estimates are meaningless — and SOL-OSI will refuse to update its beliefs from them.

121state.u = [state.u(n+1:end); dd];
122state.y = [state.y(n+1:end); dm];
126activity = min(1, (rms(state.u)/cfg.quiet_rms)^4);
141Uall = fft(state.u.*cfg.window);
142Yall = fft(state.y.*cfg.window);
145pu = abs(U).^2;
146py = abs(Y).^2;
147cyu = Y.*conj(U);
160    q = cfg.spectral_alpha;
161    state.Suu = state.Suu+q*(pu-state.Suu);
162    state.Syy = state.Syy+q*(py-state.Syy);
163    state.Syu = state.Syu+q*(cyu-state.Syu);

06 The coherent branch

How much do we trust what we hear?

With smoothed spectra in hand, the closed-loop transfer function has a classical estimator — divide the cross-spectrum by the input power, bin by bin:

(Statisticians call this the H1 estimator; it is the least-squares fit of Y on U.) The trouble is that T̂ is only as good as the data behind it. At a frequency where the drive has no energy, or where the output is mostly noise, T̂ is garbage — and one garbage bin, pushed through the min in Section 04, could fake an emergency. SOL-OSI therefore scores every bin before using it.

Score 1 — coherence: does y actually follow u?

Score 2 — leverage: was there drive there at all?

Bins more than 30 dB below the strongest bin (the power_floor) earn leverage near 0. Combining both scores with the block-level activity gate gives each bin its confidence:

Believing data exactly as much as it deserves

Here is the move that makes SOL-OSI robust rather than merely optimistic. It refuses to choose between measurement and prior knowledge — it blends them, per bin, by confidence. The prior is the healthy loop scaled by the current authority, L_prior = a L₀, pushed through Section 03’s formula to get T_prior = L_prior/(1+L_prior):

1. Confidence-weighted blend — full trust uses the measurement, zero trust falls back to the prior: T_use = confidence · T̂ + (1 − confidence) · T_prior
2. Invert Section 03 to recover the loop itself (guarding the denominator away from zero, |1 − T| ≥ 0.02): L̂ = T_use1 − T_use
3. Feed it to the chordal margin; a 3-bin moving median knocks out single-bin outliers before the min: b_st = min_keep  movmedian₃( c(L̂, −1))

This b_st is the coherent branch: the best margin estimate available when the loop is well excited. Its failure mode is equally clear — under weak drive, confidence collapses everywhere, T_use degenerates to the prior, and b_st can only repeat what it already believed. Damage that whispers instead of shouting needs a different listener.

167gamma2 = min(1, abs(state.Syu).^2./max(state.Suu.*state.Syy, eps));
168p_rel = state.Suu/max(max(state.Suu), eps);
169leverage = p_rel./(p_rel+cfg.power_floor);
170confidence = activity*gamma2.*leverage;
171That = state.Syu./max(state.Suu, eps);
172Tuse = confidence.*That+(1-confidence).*Tprior;
173den = 1-Tuse;
174small = abs(den) < 0.02;
175den(small) = 0.02*exp(1i*angle(den(small)+eps));
176Lhat = Tuse./den;
177b_st = minMasked(movmedian(chordal(Lhat, cfg.cap), 3), cfg.keep);

07 The quiet branch

Hearing damage in the noise floor

Suppose the drive is weak — the earthquake record is between pulses, activity is low, the coherent branch has gone deaf. The loop is still running, and its output still contains a faint hiss: sensor noise, ambient vibration, tiny disturbances entering everywhere. That hiss is not useless. Section 03 told us that disturbances at the output are shaped by the sensitivity — the noise floor of y carries the fingerprint of 1/(1 + L). As the loop creeps toward instability, 1 + L shrinks at the dangerous frequencies and the noise floor there bulges. The quiet branch is built to hear that bulge.

Step 1 — normalize away the healthy shape

1. Divide the measured output spectrum by the healthy sensitivity shape, precomputed from L₀: q_obs = S_yy 1/(1 + L₀)²
2. Calibrate against this specimen’s own healthy noise: the median of the first three half-second windows becomes a per-bin reference, and everything after is measured relative to it: q_obs ← q_obsq_ref ,  q_ref = median of first 3 windows

If the loop still matches its healthy self, q_obs ≈ 1 at every bin. Any systematic deviation is evidence. The branch summarizes the level in the 28–44 Hz fit band (chosen around the column resonance, where sensitivity changes show first) as the reference_ratio — you will meet it again as a safety interlock.

Step 2 — ask which damage would explain it

Now the model from the explorer earns its keep. Suppose damage acts as a loop-gain factor ρ, so the true loop is a ρ L₀. Then the noise floor should have risen by the ratio of old to new sensitivity:

SOL-OSI evaluates this on a fixed grid of 25 candidate ρ values, logarithmically spaced from 0.5 to 6, and scores each against the observation with a deliberately unglamorous loss — the median absolute log-error over the fit band:

Step 3 — believe it slowly, forget it slowly

A raw ρ_best every 0.125 s would be jumpy, so the branch keeps a damage memory: a smoothed log-domain state log ρ that moves toward the fit only in proportion to how much evidence the block carried, and leaks back toward health when there is no evidence at all:

1. Pull toward the fitted target (shape_weight ∈ [0,1] scores how decisively the evidence points at a gain change — it grows with the coherent branch’s own drive estimate and with extreme level growth): log ρ ← log ρ + activity α_res  (log ρ_target − log ρ)
2. Blind forgiveness — with no drive there is no evidence either way, so the memory decays toward ρ = 1 with a 10 s time constant: log ρ ← log ρ  (1 − α_leak(1 − activity))
3. The working damage estimate, including a fast “kick” detector explained below: ρ = e^log ρ · drive_ratio
4. And the branch’s margin: the healthy loop, re-scaled by what we currently believe about authority and damage, through the same chordal formula: b_res = min_keep  movmedian₃( c(a ρ L₀, −1))

184qobs = state.Syy./cfg.shape;
186    qobs = qobs./max(state.q_ref, eps);
187    state.reference_ratio = median(qobs(cfg.res_keep));
202qmod = abs((1+cfg.L0)./(1+a_now*cfg.rho_grid(k)*cfg.L0)).^2;
209loss(k) = median(abs(log(max(qobs(cfg.res_keep), eps))- ...
210    log(max(qmod(cfg.res_keep), eps))));
212[~, kbest] = min(loss);
216rho_drive = real(sum(weights.*conj(Lmodel).*Lhat))/ ...
217    max(sum(weights.*abs(Lmodel).^2), eps);
223raw_kick = max(rho_drive/state.rho_drive-1, 0);
224kick = min(min(raw_kick, state.ratio_prev), cfg.detector_max-1);
238state.log_rho = state.log_rho+alpha_rho*dlog_rho;
243state.log_rho = state.log_rho*(1-cfg.leak_alpha*(1-activity));
244rho = exp(state.log_rho)*drive_ratio;
245b_res = minMasked(movmedian(chordal(a_now*rho*cfg.L0, cfg.cap), 3), ...
246    cfg.keep);

08 Fusion

One number

SOL-OSI now holds two opinions: b_st, sharp but drive-hungry, and b_res, coarse but awake in the quiet. The fusion logic asks one question — who is in a position to know? — and answers it differently for driven and idle moments.

While driven: let the pessimist veto the optimist

1. How blind is the coherent branch inside the residual band right now? (drive is the slow-averaged activity) w_res = drive ·  (1 − mean_{fit band} leverage)²
2. The residual branch may only pull the margin DOWN, and only by as much as the coherent branch is blind: b_evid = b_st − w_res  max(b_st − b_res, 0)

Note the asymmetry: b_res can never raise the margin above the coherent estimate. A safety monitor accepts bad news on thin evidence and good news only on strong evidence.

While idle: report what memory justifies

With no drive and no damage memory (ρ ≈ 1), the honest idle answer is the cap — an unexcited loop shows no evidence of danger. But the gain of 32 makes the memory extraordinarily sticky: a remembered ρ only 3% above 1 already pins b_idle fully to b_res. Damage seen once is not forgotten just because the room went quiet — unless the healthy-reference interlock from Section 07 has released it (then the gain drops to 1).

Blend, smooth, clamp

248res_weight = state.drive*(1-mean(leverage(cfg.res_keep)))^2;
249b_evid = b_st-res_weight*max(b_st-b_res, 0);
250memory_gain = cfg.residual_memory_gain;
251if state.memory_released
252    memory_gain = 1;
253end
254w_mem = min(max(memory_gain*(rho-1), 0), 1);
255b_idle = cfg.cap-w_mem*(cfg.cap-b_res);
256b_raw = activity*b_evid+(1-activity)*b_idle;
257state.b = state.b+cfg.margin_alpha*(b_raw-state.b);
258b = min(max(state.b, 0), cfg.cap);

09 Closing the outer loop

Acting on it: the authority regulator

A margin nobody acts on is a seismograph in an empty room. The companion function gpt5_6_4_auth.m turns b into the knob you already used in the explorer: the authority a ∈ [0.2, 1], a gain applied to the loop. This is itself a feedback controller — a slow, deliberate outer loop wrapped around the estimator — and its target is a margin, not a position:

1. The reference margin, nudged up when the quiet branch is doing the talking (rw = residual weight from Section 08): b_ref = 0.55 + 0.02 (1 + rw)
2. The error, with a deadband so the authority does not fidget over noise-sized deviations: e_m = sign(b − b_ref) ·  max(b − b_ref − 0.008(1 + rw), 0)
3. Asymmetry: shedding authority (e < 0) acts at full strength; restoring it (e > 0) is scaled by ≈ 0.65 — fast to retreat, careful to advance: e = min(e_m, 0) + (b_ref2)²  max(e_m, 0)
4. Integrate in the log domain — so a step from 1.0 to 0.9 and from 0.4 to 0.36 are the same size move — then clamp the range and the per-update slew: log a ← log a + 5a Δt e Δa ≤ 0.05 a ∈ [0.2, 1]

And here the story closes on itself: the new a is handed back to the estimator, which uses it in the prior L_prior = a L₀, in the damage model q_mod, and in b_res. The estimator knows about the regulator’s own moves — otherwise it would diagnose its own medicine as a change in the patient.

30b_ref = prm.target+prm.guard*(1+residual_weight);
31e_raw = b-b_ref;
32deadband = prm.deadband*(1+residual_weight);
33e_margin = sign(e_raw)*max(abs(e_raw)-deadband, 0);
34recovery_weight = (b_ref*sqrt(2))^2;
35e = min(e_margin, 0)+recovery_weight*max(e_margin, 0);
36rate = prm.rate*sqrt(a0);
37log_a = log(a0)+rate*prm.dt_update*e;
38a_free = exp(min(max(log_a, log(prm.a_min)), log(prm.a_max)));
39a_next = a0+min(max(a_free-a0, -prm.max_step), prm.max_step);

10 Reference

Symbols ↔ code

Every quantity in this module, its home in sol_osi.m / gpt5_6_4_auth.m, and where it was introduced:

Provenance

The frequency-response data embedded in this page (all charts) was exported from the repository with MATLAB: C = pid(0.15, 20, 0.0015, 1/512, dt) discretized ForwardEuler at dt = 1/2048 s, multiplied by the identified plant H_dc_dm from the canonical record saved_results/tests/gpt5_6_tests/20260710_1845_…mat, sampled on cfg.freq_hz exactly as sol_osi('config', …) does. Raw export: l0_data.json beside this file. The implementation lives at control_and_HS/OSI/sol_osi.m (estimator) and control_and_HS/OSI/gpt5_6_4_auth.m (authority); the benchmark and validation notes are in research/gpt5_6_4/sol_osi/README.md.