\n');
fprintf(fid, 'Method 9.1 adaptive margin
');
fprintf(fid, 'M9.1 uses only branch 2 and branch 3: a complex closed-loop FRF where the command excites the loop, plus residual output growth for output-only evidence. All plotted and reported times are causal window-end times. Branch 1 is ignored.
\n');
writeHtmlTheory(fid, result.config);
fprintf(fid, 'True NCF references
');
fprintf(fid, '
Healthy true NCF%.3f
', ...
result.truth.healthy_ncf);
fprintf(fid, '
Damaged true NCF%.3f
', ...
result.truth.damaged_ncf);
fprintf(fid, '
Damage time%.1f s
', ...
result.truth.switch_s);
fprintf(fid, '
The dashed horizontal lines in the seed plot are these true loop-only NCF values, not seed-dependent estimates.
\n');
fprintf(fid, 'Seed spread
M9.1 seed spread includes injected noise/disturbance variability plus a bootstrap perturbation of the spectral subsegments. It is still not a full model-uncertainty interval, but it avoids plotting a deterministic averaged estimator as if finite-window spectral uncertainty were zero.
\n');
fprintf(fid, 'Between seeds, the plant, controller, damage timing, and command profile are fixed. What changes is the disturbance/noise realization, and for M9.1 the spectral-bootstrap perturbation used inside each finite-window estimate. The differences below are therefore estimator/data-window variability around the same true NCF values.
\n');
writeHtmlSeedSpread(fid, result);
fprintf(fid, 'Table comparison
\n');
writeHtmlTable(fid, result.summary);
fprintf(fid, 'Seed plots
\n', img);
fprintf(fid, 'Interpretation
M9.1 improves the three-case comparison by matching M7''s damage-reduction timing while preserving the direct phase-aware branch. The residual branch is normally input-quiet, with a frequency-matched early-warning path bounded by the direct S,T estimate; branch 1 is no longer computed or selected.
\n');
fprintf(fid, '\n\n\n');
end
function writeHtmlTheory(fid, cfg)
fprintf(fid, '\n');
fprintf(fid, '
Let \\(u=dd_{\\mathrm{lead}}\\), \\(y=dm\\), \\(d\\) be internal disturbance or noise, \\(H\\) be the plant path from command force/displacement action to measured displacement, and \\(C\\) be the controller.
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'L(\\omega,t)=C(\\omega)H(\\omega,t),\\qquad\n');
fprintf(fid, 'S(\\omega,t)=\\frac{1}{1+L(\\omega,t)},\\qquad\n');
fprintf(fid, 'T(\\omega,t)=\\frac{L(\\omega,t)}{1+L(\\omega,t)}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
The measured output is the sum of a command-correlated part and an internally excited part:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'y(\\omega,t)=T(\\omega,t)u(\\omega,t)+S(\\omega,t)d(\\omega,t)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Here \\(\\Phi_{ab}\\) means the cross-spectrum between signals \\(a\\) and \\(b\\), and \\(\\Phi_{aa}\\) is the autospectrum of \\(a\\).
\n');
fprintf(fid, '
\n');
fprintf(fid, '
Step 1: require identifiable input power
M9.1 only trusts branch 2 on bins where the input has enough power and the estimated input-output relation is coherent:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Omega_{\\mathrm{exc}}(t)=\n');
fprintf(fid, '\\left\\{\\omega:\n');
fprintf(fid, '\\Phi_{uu}(\\omega,t)\\ \\text{is large enough},\\quad\n');
fprintf(fid, '\\gamma^2(\\omega,t)\\ge \\gamma_{\\min}\n');
fprintf(fid, '\\right\\}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Step 2: estimate the command-correlated transfer
Using the output equation, the cross-spectrum with the command is:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Phi_{yu}(\\omega,t)=\n');
fprintf(fid, 'T(\\omega,t)\\Phi_{uu}(\\omega,t)\n');
fprintf(fid, '+S(\\omega,t)\\Phi_{du}(\\omega,t)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
If the internal disturbance is not correlated with the command, \\(\\Phi_{du}\\approx0\\). Then:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\widehat T(\\omega,t)=\n');
fprintf(fid, '\\frac{\\Phi_{yu}(\\omega,t)}{\\Phi_{uu}(\\omega,t)}\n');
fprintf(fid, '\\approx T(\\omega,t)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Step 3: convert \\(S,T\\) into margin
For unity feedback, \\(S=1-T\\). The NCF loop distance at one frequency can be written directly from \\(S,T\\):
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'b(\\omega,t)=\n');
fprintf(fid, '\\frac{1}{\\sqrt{2}\\sqrt{|S(\\omega,t)|^2+|T(\\omega,t)|^2}}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
M9.1 uses the largest distance over qualified excited/coherent bins:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\tilde b_{ST}(t)=\n');
fprintf(fid, '\\frac{1}{\\sqrt{2}\\,\n');
fprintf(fid, '\\max_{\\omega\\in\\Omega_{\\mathrm{exc}}(t)}\n');
fprintf(fid, '\\left(\\sqrt{|\\widehat S(\\omega,t)|^2+|\\widehat T(\\omega,t)|^2}\\right)}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
A small finite-window correction sharpens sub-nominal branch-2 margins while leaving the healthy nominal point anchored:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'b_{ST}(t)=\n');
fprintf(fid, 'b_0\\left(\\frac{\\tilde b_{ST}(t)}{b_0}\\right)^{%.2f},\\qquad \\tilde b_{ST}(t)
\n');
fprintf(fid, '
This branch is the phase-aware branch. If the dangerous frequency is externally excited, it can detect margin loss caused by phase delay even before a large output-only resonance is visible.
\n');
fprintf(fid, '
\n');
fprintf(fid, '
Branch 3: residual output evidence
\n
\n');
fprintf(fid, '
Step 1: decompose the output spectrum
When \\(u\\) and \\(d\\) are approximately uncorrelated, the total output spectrum is:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Phi_{yy}(\\omega,t)=\n');
fprintf(fid, '\\underbrace{|T(\\omega,t)|^2\\Phi_{uu}(\\omega,t)}_{\\text{command-correlated output}}\n');
fprintf(fid, '+\n');
fprintf(fid, '\\underbrace{|S(\\omega,t)|^2\\Phi_{dd}(\\omega)}_{\\text{residual/internal output}}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Step 2: use coherence to remove the command-correlated part
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\gamma^2(\\omega,t)=\n');
fprintf(fid, '\\frac{|\\Phi_{yu}(\\omega,t)|^2}\n');
fprintf(fid, '{\\Phi_{uu}(\\omega,t)\\Phi_{yy}(\\omega,t)}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
The residual output spectrum is approximated as:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Phi_{\\mathrm{res}}(\\omega,t)\\approx\n');
fprintf(fid, '\\Phi_{yy}(\\omega,t)\\left(1-\\gamma^2(\\omega,t)\\right)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Step 3: search residual growth
Branch 3 compares the current residual spectrum with the healthy reference spectrum. The normal fallback searches outside the branch-2 excited/coherent bins:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Omega_{\\mathrm{quiet}}(t)=\\Omega_{\\mathrm{wide}}\\setminus\\Omega_{\\mathrm{exc}}(t),\\qquad\n');
fprintf(fid, 'q(t)=\n');
fprintf(fid, '\\max_{\\omega\\in\\Omega_{\\mathrm{quiet}}(t)}\n');
fprintf(fid, '\\frac{\\Phi_{\\mathrm{res}}(\\omega,t)}\n');
fprintf(fid, '{\\Phi_{\\mathrm{res,ref}}(\\omega)}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
Then it maps the residual growth to a conservative margin proxy:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'b_{\\mathrm{wide}}(t)=\n');
fprintf(fid, 'b_0(\\omega^*)\\sqrt{\\frac{q_{\\mathrm{ref}}}{q(t)}}\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
For faster damage-reduction detection, M9.1 also checks a full-band residual candidate. It is accepted only when it binds within %.1f Hz of the direct S,T branch and S,T has already moved below nominal but not below the relative alarm level. In that case the residual candidate is bounded below by %.0f%% of the direct S,T margin.
\n', cfg.m9_res_early_match_hz, 100*cfg.m9_res_early_st_floor_frac);
fprintf(fid, '
Why this sees what branch 2 can miss
If a mode grows from internal disturbance or noise, then the residual spectrum contains:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Phi_{\\mathrm{res}}(\\omega,t)\\approx\n');
fprintf(fid, '|S(\\omega,t)|^2\\Phi_{dd}(\\omega)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
So branch 3 can catch output-only resonance growth when the command has almost no power at that frequency, and the bounded early-warning path can also catch the first residual growth at an already-observed S,T frequency:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, '\\Phi_{uu}(\\omega)\\approx0\n');
fprintf(fid, '\\]
\n
\n');
fprintf(fid, '
Final branch selection
\n
\n');
fprintf(fid, '
M9.1 computes no known-mode peak/floor branch. Each window keeps only the lower margin from branch 2 and branch 3:
\n');
fprintf(fid, '
\\[\n');
fprintf(fid, 'b_{M9.1}(t)=\\min\\left(b_{ST}(t), b_{\\mathrm{wide}}(t)\\right)\n');
fprintf(fid, '\\]
\n');
fprintf(fid, '
This makes branch 3 the only output-only fallback. There is no branch-1 peak mask, no branch-1 calibrated proxy, and no final minimum against a branch-1 estimate.
\n');
fprintf(fid, '
\n');
end
function writeSummaryTable(fid, summary)
fprintf(fid, '| Case | M7 healthy | M7 damaged | M7 abs error | M7 seeds | M7 median s | M9.1 healthy | M9.1 damaged | M9.1 abs error | M9.1 seeds | M9.1 median s | M9.1 freq Hz | M7 ms | M9.1 ms | M9.1/M7 |\n');
fprintf(fid, '|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|\n');
for i = 1:height(summary)
fprintf(fid, '| %s | %.3f | %.3f | %.3f | %d | %.1f | %.3f | %.3f | %.3f | %d | %.1f | %.1f | %.0f | %.0f | %.2f |\n', ...
summary.Case(i), summary.M7Healthy(i), summary.M7Damaged(i), ...
summary.M7DamageAbsError(i), summary.M7DetectSeeds(i), summary.M7MedianDetect_s(i), ...
summary.M9Healthy(i), summary.M9Damaged(i), summary.M9DamageAbsError(i), ...
summary.M9DetectSeeds(i), summary.M9MedianDetect_s(i), summary.M9MedianFreqHz(i), ...
summary.M7RuntimeMs(i), summary.M9RuntimeMs(i), summary.M9RuntimeRatio(i));
end
end
function writeHtmlTable(fid, summary)
fprintf(fid, '
| Case | M7 damaged | M7 err | M7 seeds | M9.1 damaged | M9.1 err | M9.1 seeds | M9.1 freq | M7 ms | M9.1 ms | M9.1/M7 |
|---|
\n');
for i = 1:height(summary)
cls = "good";
if summary.M9DamageAbsError(i) > summary.M7DamageAbsError(i)
cls = "bad";
end
speed_cls = "good";
if summary.M9RuntimeRatio(i) > 1.20
speed_cls = "bad";
end
fprintf(fid, '| %s | %.3f | %.3f | %d/10 | %.3f | %.3f | %d/10 | %.1f Hz | %.0f | %.0f | %.2fx |
\n', ...
summary.Case(i), summary.M7Damaged(i), summary.M7DamageAbsError(i), ...
summary.M7DetectSeeds(i), cls, summary.M9Damaged(i), cls, ...
summary.M9DamageAbsError(i), cls, summary.M9DetectSeeds(i), summary.M9MedianFreqHz(i), ...
summary.M7RuntimeMs(i), speed_cls, summary.M9RuntimeMs(i), speed_cls, summary.M9RuntimeRatio(i));
end
fprintf(fid, '
\n');
end
function writeHtmlSeedSpread(fid, result)
fprintf(fid, '
| Case | What differs | M9.1 damaged seed range | M9.1 max seed spread | M9.1 detection range | M7 damaged seed range | M7 detection range |
|---|
\n');
for i = 1:numel(result.cases)
c = result.cases(i);
d9 = c.tc9 >= result.config.switch_s + result.config.m9_buf;
d7 = c.tc7 >= result.config.switch_s + result.config.m7_buf;
m9_seed_damaged = mean(c.B9(d9, :), 1, 'omitnan');
m7_seed_damaged = mean(c.B7(d7, :), 1, 'omitnan');
m9_window_spread = max(c.B9, [], 2, 'omitnan') ...
- min(c.B9, [], 2, 'omitnan');
m9_detect = formatRange(c.td9);
m7_detect = formatRange(c.td7);
fprintf(fid, '| %s | noise realization + M9.1 spectral bootstrap | %.3f-%.3f | %.3f | %s | %.3f-%.3f | %s |
\n', ...
c.name, min(m9_seed_damaged, [], 'omitnan'), ...
max(m9_seed_damaged, [], 'omitnan'), ...
max(m9_window_spread, [], 'omitnan'), m9_detect, ...
min(m7_seed_damaged, [], 'omitnan'), ...
max(m7_seed_damaged, [], 'omitnan'), m7_detect);
end
fprintf(fid, '
\n');
fprintf(fid, '
Damaged seed range is the per-seed mean after the estimator buffer has cleared the damage step. Max seed spread is the largest vertical separation between M9.1 seed traces at any plotted timestamp.
\n');
end
function text = formatRange(values)
values = values(~isnan(values));
if isempty(values)
text = 'none';
elseif min(values) == max(values)
text = sprintf('%.1f s', min(values));
else
text = sprintf('%.1f-%.1f s', min(values), max(values));
end
end
function dd = grab(file, t0, N)
R = load(file);
i0 = find(R.result.t >= t0, 1);
dd = R.result.dd(i0:i0+N-1);
dd = dd - mean(dd);
end
function b = minNcf(C, H, band)
f = logspace(log10(band(1)), log10(band(2)), 4000).';
L = squeeze(freqresp(C*H, 2*pi*f));
b = min(abs(1 + L)./sqrt(2*(1 + abs(L).^2)));
end
function td = relhit(ncf, tc, thr, consec)
below = ncf < thr;
below(isnan(ncf)) = false;
hit = find(movsum(below, [consec-1, 0]) >= consec, 1);
td = NaN;
if ~isempty(hit)
td = tc(hit);
end
end