Example: Multi-Trajectory Ensemble Averaging¶

Demonstrates how multiple independent trajectories improve spectral estimation by ensemble averaging. The same physical plant is measured L times with different excitation and noise realizations; averaging reduces variance by a factor of L without sacrificing frequency resolution.

Functions demonstrated: sid.freq_bt, sid.freq_map, sid.spectrogram, sid.ltv_disc, sid.ltv_disc_frozen.

Plants. Sections 1, 2, and 4 use a 2-mass spring-damper chain (m=[1,1], k=[100,80], c=[2,2]). Section 3 uses the SDOF plant from example_spectrogram (m=1, k=4·10⁴, c=20) driven by a chirp force buried in noise.

Translated from exampleMultiTrajectory.m.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd, util_msd_ltv
import sid

%matplotlib inline

import numpy as np import matplotlib.pyplot as plt from util_msd import util_msd, util_msd_ltv import sid %matplotlib inline

1. LTI ensemble averaging — tighter confidence bands¶

Generate L = 10 trajectories of the 2-mass SMD. Compare sid.freq_bt confidence bands using 1 vs 10 trajectories.

In [2]:

rng = np.random.default_rng(5001)

m  = np.array([1.0, 1.0])
k  = np.array([100.0, 80.0])
c  = np.array([2.0, 2.0])
F  = np.array([[1.0], [0.0]])
Ts = 0.01
N, L = 2000, 10

Ad, Bd = util_msd(m, k, c, F, Ts)

y_all = np.zeros((N, 1, L))
u_all = np.zeros((N, 1, L))
for l in range(L):
    u_all[:, 0, l] = rng.standard_normal(N)
    xs = np.zeros((N + 1, 4))
    for step in range(N):
        xs[step + 1] = Ad @ xs[step] + Bd[:, 0] * u_all[step, 0, l]
    y_all[:, 0, l] = xs[1:, 0] + 5e-4 * rng.standard_normal(N)

# Single trajectory
r1 = sid.freq_bt(y_all[:, :, 0], u_all[:, :, 0], window_size=80, sample_time=Ts)

# Multi-trajectory ensemble
rL = sid.freq_bt(y_all, u_all, window_size=80, sample_time=Ts)

print(f'Single trajectory: max response_std = {r1.response_std.max():.3e}')
print(f'{L} trajectories:   max response_std = {rL.response_std.max():.3e}')
print(f'Ratio: {rL.response_std.max() / r1.response_std.max():.2f} '
      f'(expected ~{1 / np.sqrt(L):.2f} = 1/sqrt({L}))')

fig, axes = plt.subplots(2, 1, figsize=(8, 8))
sid.bode_plot(r1, confidence=3, ax=(axes[0], axes[0].twinx()))
axes[0].set_title('Single trajectory')
sid.bode_plot(rL, confidence=3, ax=(axes[1], axes[1].twinx()))
axes[1].set_title(f'Ensemble of {L} trajectories')
plt.tight_layout()
plt.show()

rng = np.random.default_rng(5001) m = np.array([1.0, 1.0]) k = np.array([100.0, 80.0]) c = np.array([2.0, 2.0]) F = np.array([[1.0], [0.0]]) Ts = 0.01 N, L = 2000, 10 Ad, Bd = util_msd(m, k, c, F, Ts) y_all = np.zeros((N, 1, L)) u_all = np.zeros((N, 1, L)) for l in range(L): u_all[:, 0, l] = rng.standard_normal(N) xs = np.zeros((N + 1, 4)) for step in range(N): xs[step + 1] = Ad @ xs[step] + Bd[:, 0] * u_all[step, 0, l] y_all[:, 0, l] = xs[1:, 0] + 5e-4 * rng.standard_normal(N) # Single trajectory r1 = sid.freq_bt(y_all[:, :, 0], u_all[:, :, 0], window_size=80, sample_time=Ts) # Multi-trajectory ensemble rL = sid.freq_bt(y_all, u_all, window_size=80, sample_time=Ts) print(f'Single trajectory: max response_std = {r1.response_std.max():.3e}') print(f'{L} trajectories: max response_std = {rL.response_std.max():.3e}') print(f'Ratio: {rL.response_std.max() / r1.response_std.max():.2f} ' f'(expected ~{1 / np.sqrt(L):.2f} = 1/sqrt({L}))') fig, axes = plt.subplots(2, 1, figsize=(8, 8)) sid.bode_plot(r1, confidence=3, ax=(axes[0], axes[0].twinx())) axes[0].set_title('Single trajectory') sid.bode_plot(rL, confidence=3, ax=(axes[1], axes[1].twinx())) axes[1].set_title(f'Ensemble of {L} trajectories') plt.tight_layout() plt.show()

Single trajectory: max response_std = 2.839e-03
10 trajectories:   max response_std = 8.941e-04
Ratio: 0.31 (expected ~0.32 = 1/sqrt(10))

2. LTV time-varying map — sharper transition detection¶

Same 2-mass plant but with k₁ snapping from 200 N/m to 50 N/m at the mid-point of the record — a sudden structural change. The time-frequency map should show the transition; with more trajectories, the map gets crisper because segment-level variance is averaged away.

In [3]:

rng2 = np.random.default_rng(5002)

N_tv, L_tv = 4000, 5
k_step = np.zeros((2, N_tv))
k_step[0, :N_tv // 2] = 200.0
k_step[0, N_tv // 2:] = 50.0
k_step[1, :] = 80.0
m_tv = np.broadcast_to(m[:, None], (2, N_tv))
c_tv = np.broadcast_to(c[:, None], (2, N_tv))

Ad_tv, Bd_tv = util_msd_ltv(m_tv, k_step, c_tv, F, Ts)

y_ltv = np.zeros((N_tv, 1, L_tv))
u_ltv = np.zeros((N_tv, 1, L_tv))
for l in range(L_tv):
    u_ltv[:, 0, l] = rng2.standard_normal(N_tv)
    xs = np.zeros((N_tv + 1, 4))
    for step in range(N_tv):
        xs[step + 1] = (Ad_tv[:, :, step] @ xs[step]
                        + Bd_tv[:, 0, step] * u_ltv[step, 0, l])
    y_ltv[:, 0, l] = xs[1:, 0] + 5e-4 * rng2.standard_normal(N_tv)

r1_map = sid.freq_map(y_ltv[:, :, 0], u_ltv[:, :, 0], segment_length=256, sample_time=Ts)
rL_map = sid.freq_map(y_ltv, u_ltv, segment_length=256, sample_time=Ts)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.map_plot(r1_map, plot_type='magnitude', ax=axes[0])
axes[0].set_title('Single trajectory')
sid.map_plot(rL_map, plot_type='magnitude', ax=axes[1])
axes[1].set_title(f'Ensemble of {L_tv} trajectories')
plt.tight_layout()
plt.show()

rng2 = np.random.default_rng(5002) N_tv, L_tv = 4000, 5 k_step = np.zeros((2, N_tv)) k_step[0, :N_tv // 2] = 200.0 k_step[0, N_tv // 2:] = 50.0 k_step[1, :] = 80.0 m_tv = np.broadcast_to(m[:, None], (2, N_tv)) c_tv = np.broadcast_to(c[:, None], (2, N_tv)) Ad_tv, Bd_tv = util_msd_ltv(m_tv, k_step, c_tv, F, Ts) y_ltv = np.zeros((N_tv, 1, L_tv)) u_ltv = np.zeros((N_tv, 1, L_tv)) for l in range(L_tv): u_ltv[:, 0, l] = rng2.standard_normal(N_tv) xs = np.zeros((N_tv + 1, 4)) for step in range(N_tv): xs[step + 1] = (Ad_tv[:, :, step] @ xs[step] + Bd_tv[:, 0, step] * u_ltv[step, 0, l]) y_ltv[:, 0, l] = xs[1:, 0] + 5e-4 * rng2.standard_normal(N_tv) r1_map = sid.freq_map(y_ltv[:, :, 0], u_ltv[:, :, 0], segment_length=256, sample_time=Ts) rL_map = sid.freq_map(y_ltv, u_ltv, segment_length=256, sample_time=Ts) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) sid.map_plot(r1_map, plot_type='magnitude', ax=axes[0]) axes[0].set_title('Single trajectory') sid.map_plot(rL_map, plot_type='magnitude', ax=axes[1]) axes[1].set_title(f'Ensemble of {L_tv} trajectories') plt.tight_layout() plt.show()

3. Spectrogram averaging — chirp in noise¶

A chirp force driving an SDOF plant, with heavy measurement noise added to the response. Single trajectory: chirp barely visible above the noise floor. L = 8 trajectories averaged: the chirp track emerges clearly.

In [4]:

rng3 = np.random.default_rng(5003)

N_sp, L_sp, Fs = 2000, 8, 1000
Ts_sp = 1 / Fs
t_sp = np.arange(N_sp) * Ts_sp

# SDOF plant with resonance at 31.83 Hz (same as example_spectrogram)
m_sp  = np.array([1.0])
k_sp  = np.array([4e4])
c_sp  = np.array([20.0])
F_sp  = np.array([[1.0]])
Ad_sp, Bd_sp = util_msd(m_sp, k_sp, c_sp, F_sp, Ts_sp)

# Chirp force: 20 -> 60 Hz
f0, f1 = 20, 60
u_chirp = np.cos(2 * np.pi * (f0 + (f1 - f0) / (2 * t_sp[-1]) * t_sp) * t_sp)

x_all = np.zeros((N_sp, 1, L_sp))
for l in range(L_sp):
    xs = np.zeros((N_sp + 1, 2))
    for step in range(N_sp):
        xs[step + 1] = Ad_sp @ xs[step] + Bd_sp[:, 0] * u_chirp[step]
    # Add noise that is LARGE relative to the tiny response amplitude
    x_all[:, 0, l] = xs[1:, 0] + 1e-4 * rng3.standard_normal(N_sp)

r1_spec = sid.spectrogram(x_all[:, :, 0], window_length=128, sample_time=Ts_sp)
rL_spec = sid.spectrogram(x_all,          window_length=128, sample_time=Ts_sp)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.spectrogram_plot(r1_spec, ax=axes[0])
axes[0].set_title('Single trajectory')
sid.spectrogram_plot(rL_spec, ax=axes[1])
axes[1].set_title(f'Ensemble of {L_sp} trajectories')
plt.tight_layout()
plt.show()

rng3 = np.random.default_rng(5003) N_sp, L_sp, Fs = 2000, 8, 1000 Ts_sp = 1 / Fs t_sp = np.arange(N_sp) * Ts_sp # SDOF plant with resonance at 31.83 Hz (same as example_spectrogram) m_sp = np.array([1.0]) k_sp = np.array([4e4]) c_sp = np.array([20.0]) F_sp = np.array([[1.0]]) Ad_sp, Bd_sp = util_msd(m_sp, k_sp, c_sp, F_sp, Ts_sp) # Chirp force: 20 -> 60 Hz f0, f1 = 20, 60 u_chirp = np.cos(2 * np.pi * (f0 + (f1 - f0) / (2 * t_sp[-1]) * t_sp) * t_sp) x_all = np.zeros((N_sp, 1, L_sp)) for l in range(L_sp): xs = np.zeros((N_sp + 1, 2)) for step in range(N_sp): xs[step + 1] = Ad_sp @ xs[step] + Bd_sp[:, 0] * u_chirp[step] # Add noise that is LARGE relative to the tiny response amplitude x_all[:, 0, l] = xs[1:, 0] + 1e-4 * rng3.standard_normal(N_sp) r1_spec = sid.spectrogram(x_all[:, :, 0], window_length=128, sample_time=Ts_sp) rL_spec = sid.spectrogram(x_all, window_length=128, sample_time=Ts_sp) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) sid.spectrogram_plot(r1_spec, ax=axes[0]) axes[0].set_title('Single trajectory') sid.spectrogram_plot(rL_spec, ax=axes[1]) axes[1].set_title(f'Ensemble of {L_sp} trajectories') plt.tight_layout() plt.show()

4. COSMIC + sid.freq_map consistency¶

Use the same multi-trajectory dataset to run both COSMIC (state-space) and sid.freq_map (non-parametric). The frozen transfer function from COSMIC should match the ensemble-averaged spectral map, since they are estimating the same underlying plant from the same data.

For this test we go back to the 1-DoF SMD with ramping stiffness from example_ltv_disc, but use fewer trajectories and a shorter record to keep runtime low.

In [5]:

rng4 = np.random.default_rng(5004)

m_s  = np.array([1.0])
c_s  = np.array([2.0])
F_s  = np.array([[1.0]])
p, q = 2, 1
N_co, L_co = 80, 10

k_co = np.linspace(200.0, 50.0, N_co).reshape(1, N_co)
m_co = np.broadcast_to(m_s[:, None], (1, N_co))
c_co = np.broadcast_to(c_s[:, None], (1, N_co))
Ad_co, Bd_co = util_msd_ltv(m_co, k_co, c_co, F_s, Ts)

X = np.zeros((N_co + 1, p, L_co))
U = rng4.standard_normal((N_co, q, L_co))
for l in range(L_co):
    X[0, :, l] = rng4.standard_normal(p)
    for step in range(N_co):
        X[step + 1, :, l] = (Ad_co[:, :, step] @ X[step, :, l]
                             + Bd_co[:, :, step] @ U[step, :, l]
                             + 0.01 * rng4.standard_normal(p))

ltv = sid.ltv_disc(X, U, lambda_='auto', uncertainty=True)
print(f'COSMIC identified A(0) =\n{ltv.a[:, :, 0]}')
print(f'COSMIC identified A(N-1) =\n{ltv.a[:, :, -1]}')

# Ensemble-averaged freq_map using position (x1) as output
y_freq = X[:N_co, 0:1, :].copy()
u_freq = U.copy()
fmap = sid.freq_map(y_freq, u_freq, segment_length=min(N_co, 30), sample_time=Ts)

print(f'\nfreq_map ensemble: {fmap.num_trajectories} trajectories, '
      f'{len(fmap.time)} time segments.')
print(f'Both COSMIC and freq_map use the same {L_co}-trajectory dataset.')

rng4 = np.random.default_rng(5004) m_s = np.array([1.0]) c_s = np.array([2.0]) F_s = np.array([[1.0]]) p, q = 2, 1 N_co, L_co = 80, 10 k_co = np.linspace(200.0, 50.0, N_co).reshape(1, N_co) m_co = np.broadcast_to(m_s[:, None], (1, N_co)) c_co = np.broadcast_to(c_s[:, None], (1, N_co)) Ad_co, Bd_co = util_msd_ltv(m_co, k_co, c_co, F_s, Ts) X = np.zeros((N_co + 1, p, L_co)) U = rng4.standard_normal((N_co, q, L_co)) for l in range(L_co): X[0, :, l] = rng4.standard_normal(p) for step in range(N_co): X[step + 1, :, l] = (Ad_co[:, :, step] @ X[step, :, l] + Bd_co[:, :, step] @ U[step, :, l] + 0.01 * rng4.standard_normal(p)) ltv = sid.ltv_disc(X, U, lambda_='auto', uncertainty=True) print(f'COSMIC identified A(0) =\n{ltv.a[:, :, 0]}') print(f'COSMIC identified A(N-1) =\n{ltv.a[:, :, -1]}') # Ensemble-averaged freq_map using position (x1) as output y_freq = X[:N_co, 0:1, :].copy() u_freq = U.copy() fmap = sid.freq_map(y_freq, u_freq, segment_length=min(N_co, 30), sample_time=Ts) print(f'\nfreq_map ensemble: {fmap.num_trajectories} trajectories, ' f'{len(fmap.time)} time segments.') print(f'Both COSMIC and freq_map use the same {L_co}-trajectory dataset.')

COSMIC identified A(0) =
[[ 0.98833003  0.00746722]
 [-1.96185775  0.97388733]]
COSMIC identified A(N-1) =
[[ 0.98366657  0.0128609 ]
 [-0.64828547  0.99797256]]

freq_map ensemble: 10 trajectories, 4 time segments.
Both COSMIC and freq_map use the same 10-trajectory dataset.