Example: LTV State-Space Identification with `sid.ltv_disc`¶

This example demonstrates sid.ltv_disc (the COSMIC algorithm) for identifying time-varying discrete linear systems of the form

$$x(k+1) = A(k)\,x(k) + B(k)\,u(k)$$

from multiple trajectories. We also exercise sid.ltv_disc_tune (regularisation tuning) and sid.ltv_disc_frozen (frozen frequency response at selected time steps).

Plant. A 1-DoF spring-mass-damper, written as a state-space model with state [position; velocity]:

Parameter	Value
`m`	`1 kg`
`k` (baseline)	`100 N/m`
`c`	`2 N·s/m`
`Ts`	`0.01 s`

With these parameters the discrete-time dynamics matrix at the baseline stiffness is

Ad ≈ [[ 0.995,  0.010],
      [-0.988,  0.975]]

The A[1, 0] entry is approximately -k·Ts (force-to-velocity coupling), so when we ramp k(t) later, the recovered A[1, 0](k) should track -k(t)·Ts.

At the end we exercise COSMIC on a nonlinear Duffing oscillator (util_msd_nl) and recover the amplitude-dependent local linearisation — the classical "weakly nonlinear → apparently LTV" interpretation of COSMIC.

Translated from exampleLTVdisc.m.

In [1]:

Copied!





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

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

%matplotlib inline

1. LTI system recovery¶

First, verify that sid.ltv_disc correctly identifies a known LTI system by passing multiple trajectories of the constant-stiffness 1-DoF plant. The mean of the recovered A(k) across time should match the true Ad to within the process-noise-limited resolution.

In [2]:

Copied!





rng = np.random.default_rng(100)

m  = np.array([1.0])
k  = np.array([100.0])
c  = np.array([2.0])
F  = np.array([[1.0]])
Ts = 0.01

Ad_true, Bd_true = util_msd(m, k, c, F, Ts)   # 2x2, 2x1
p, q = 2, 1

N, L = 50, 10
sigma = 0.01

X = np.zeros((N + 1, p, L))
U = rng.standard_normal((N, q, L))
for l in range(L):
    X[0, :, l] = rng.standard_normal(p)
    for step in range(N):
        X[step + 1, :, l] = (Ad_true @ X[step, :, l]
                             + Bd_true @ U[step, :, l]
                             + sigma * rng.standard_normal(p))

result_lti = sid.ltv_disc(X, U, lambda_=1e5)
A_mean = np.mean(result_lti.a, axis=2)
print('True Ad:')
print(Ad_true)
print('Mean recovered A(k):')
print(A_mean)
print(f'Recovery error (Frobenius): {np.linalg.norm(A_mean - Ad_true):.4e}')
rng = np.random.default_rng(100)

m  = np.array([1.0])
k  = np.array([100.0])
c  = np.array([2.0])
F  = np.array([[1.0]])
Ts = 0.01

Ad_true, Bd_true = util_msd(m, k, c, F, Ts)   # 2x2, 2x1
p, q = 2, 1

N, L = 50, 10
sigma = 0.01

X = np.zeros((N + 1, p, L))
U = rng.standard_normal((N, q, L))
for l in range(L):
    X[0, :, l] = rng.standard_normal(p)
    for step in range(N):
        X[step + 1, :, l] = (Ad_true @ X[step, :, l]
                             + Bd_true @ U[step, :, l]
                             + sigma * rng.standard_normal(p))

result_lti = sid.ltv_disc(X, U, lambda_=1e5)
A_mean = np.mean(result_lti.a, axis=2)
print('True Ad:')
print(Ad_true)
print('Mean recovered A(k):')
print(A_mean)
print(f'Recovery error (Frobenius): {np.linalg.norm(A_mean - Ad_true):.4e}')

True Ad:
[[ 0.9950373   0.00988417]
 [-0.98841706  0.97526896]]
Mean recovered A(k):
[[ 0.99419405  0.00983385]
 [-0.98881978  0.97516648]]
Recovery error (Frobenius): 9.4143e-04

2. LTV system: time-varying stiffness¶

Now the spring constant ramps linearly from k₁ = 200 N/m to k₁ = 50 N/m over the record. We build the per-step discrete dynamics stack with util_msd_ltv and simulate the LTV recursion. Because Ad[1, 0] ≈ -k·Ts for small Ts, the recovered A[1, 0](k) should drift from about -2.0 to -0.5.

In [3]:

Copied!





rng2 = np.random.default_rng(200)

N_tv, L_tv = 250, 30
k_tv = np.linspace(200.0, 50.0, N_tv).reshape(1, N_tv)
m_tv = np.broadcast_to(m[:, None], (1, N_tv))
c_tv = np.broadcast_to(c[:, None], (1, N_tv))
Ad_tv, Bd_tv = util_msd_ltv(m_tv, k_tv, c_tv, F, Ts)   # (2,2,N), (2,1,N)

X_tv = np.zeros((N_tv + 1, p, L_tv))
U_tv = rng2.standard_normal((N_tv, q, L_tv))
for l in range(L_tv):
    X_tv[0, :, l] = rng2.standard_normal(p)
    for step in range(N_tv):
        X_tv[step + 1, :, l] = (Ad_tv[:, :, step] @ X_tv[step, :, l]
                                + Bd_tv[:, :, step] @ U_tv[step, :, l]
                                + 0.01 * rng2.standard_normal(p))

result_tv = sid.ltv_disc(X_tv, U_tv, lambda_='auto')
print(f'auto lambda: {result_tv.lambda_[0]:.3e}')

fig, ax = plt.subplots(figsize=(8, 5))
kk = np.arange(N_tv)
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True A[1, 0](k) ≈ -k(t)·Ts')
ax.plot(kk, result_tv.a[1, 0, :], 'b', label='COSMIC (auto λ)')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('LTV identification: recovered A[1, 0] tracks the ramping stiffness')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()
rng2 = np.random.default_rng(200)

N_tv, L_tv = 250, 30
k_tv = np.linspace(200.0, 50.0, N_tv).reshape(1, N_tv)
m_tv = np.broadcast_to(m[:, None], (1, N_tv))
c_tv = np.broadcast_to(c[:, None], (1, N_tv))
Ad_tv, Bd_tv = util_msd_ltv(m_tv, k_tv, c_tv, F, Ts)   # (2,2,N), (2,1,N)

X_tv = np.zeros((N_tv + 1, p, L_tv))
U_tv = rng2.standard_normal((N_tv, q, L_tv))
for l in range(L_tv):
    X_tv[0, :, l] = rng2.standard_normal(p)
    for step in range(N_tv):
        X_tv[step + 1, :, l] = (Ad_tv[:, :, step] @ X_tv[step, :, l]
                                + Bd_tv[:, :, step] @ U_tv[step, :, l]
                                + 0.01 * rng2.standard_normal(p))

result_tv = sid.ltv_disc(X_tv, U_tv, lambda_='auto')
print(f'auto lambda: {result_tv.lambda_[0]:.3e}')

fig, ax = plt.subplots(figsize=(8, 5))
kk = np.arange(N_tv)
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True A[1, 0](k) ≈ -k(t)·Ts')
ax.plot(kk, result_tv.a[1, 0, :], 'b', label='COSMIC (auto λ)')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('LTV identification: recovered A[1, 0] tracks the ramping stiffness')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

auto lambda: 3.728e-01

No description has been provided for this image

3. Multi-trajectory benefit¶

More trajectories provide more information and reduce estimation error. Compare L = 3 trajectories against L = 20.

In [4]:

Copied!





rng3 = np.random.default_rng(300)
L_few, L_many = 5, 30

X_few = X_tv[:, :, :L_few]
U_few = U_tv[:, :, :L_few]

X_many = np.zeros((N_tv + 1, p, L_many))
U_many = rng3.standard_normal((N_tv, q, L_many))
for l in range(L_many):
    X_many[0, :, l] = rng3.standard_normal(p)
    for step in range(N_tv):
        X_many[step + 1, :, l] = (Ad_tv[:, :, step] @ X_many[step, :, l]
                                  + Bd_tv[:, :, step] @ U_many[step, :, l]
                                  + 0.01 * rng3.standard_normal(p))

r_few  = sid.ltv_disc(X_few,  U_few,  lambda_='auto')
r_many = sid.ltv_disc(X_many, U_many, lambda_='auto')

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True')
ax.plot(kk, r_few.a[1, 0, :],  'r', label=f'L = {L_few}')
ax.plot(kk, r_many.a[1, 0, :], 'b', label=f'L = {L_many}')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('Effect of number of trajectories')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()
rng3 = np.random.default_rng(300)
L_few, L_many = 5, 30

X_few = X_tv[:, :, :L_few]
U_few = U_tv[:, :, :L_few]

X_many = np.zeros((N_tv + 1, p, L_many))
U_many = rng3.standard_normal((N_tv, q, L_many))
for l in range(L_many):
    X_many[0, :, l] = rng3.standard_normal(p)
    for step in range(N_tv):
        X_many[step + 1, :, l] = (Ad_tv[:, :, step] @ X_many[step, :, l]
                                  + Bd_tv[:, :, step] @ U_many[step, :, l]
                                  + 0.01 * rng3.standard_normal(p))

r_few  = sid.ltv_disc(X_few,  U_few,  lambda_='auto')
r_many = sid.ltv_disc(X_many, U_many, lambda_='auto')

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True')
ax.plot(kk, r_few.a[1, 0, :],  'r', label=f'L = {L_few}')
ax.plot(kk, r_many.a[1, 0, :], 'b', label=f'L = {L_many}')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('Effect of number of trajectories')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

4. Validation-based lambda tuning with `sid.ltv_disc_tune`¶

Split trajectories into training and validation sets, then search over a grid of lambda values to minimise validation prediction error.

In [5]:

Copied!





rng4 = np.random.default_rng(400)
L_total, L_train = 30, 20

X_all = np.zeros((N_tv + 1, p, L_total))
U_all = rng4.standard_normal((N_tv, q, L_total))
for l in range(L_total):
    X_all[0, :, l] = rng4.standard_normal(p)
    for step in range(N_tv):
        X_all[step + 1, :, l] = (Ad_tv[:, :, step] @ X_all[step, :, l]
                                 + Bd_tv[:, :, step] @ U_all[step, :, l]
                                 + 0.01 * rng4.standard_normal(p))

X_train = X_all[:, :, :L_train]
U_train = U_all[:, :, :L_train]
X_val   = X_all[:, :, L_train:]
U_val   = U_all[:, :, L_train:]

grid_lam = np.logspace(-3, 6, 30)
best_result, best_lambda, all_losses = sid.ltv_disc_tune(
    X_train, U_train, X_val, U_val,
    lambda_grid=grid_lam,
)
print(f'Validation-tuned lambda: {best_lambda:.2e}')

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(grid_lam, all_losses, 'b.-')
ax.semilogx(best_lambda, np.min(all_losses), 'ro', ms=10, lw=2)
ax.set_xlabel(r'$\lambda$')
ax.set_ylabel('Validation RMSE')
ax.set_title('Lambda tuning: validation loss curve')
ax.grid(True)
plt.tight_layout()
plt.show()
rng4 = np.random.default_rng(400)
L_total, L_train = 30, 20

X_all = np.zeros((N_tv + 1, p, L_total))
U_all = rng4.standard_normal((N_tv, q, L_total))
for l in range(L_total):
    X_all[0, :, l] = rng4.standard_normal(p)
    for step in range(N_tv):
        X_all[step + 1, :, l] = (Ad_tv[:, :, step] @ X_all[step, :, l]
                                 + Bd_tv[:, :, step] @ U_all[step, :, l]
                                 + 0.01 * rng4.standard_normal(p))

X_train = X_all[:, :, :L_train]
U_train = U_all[:, :, :L_train]
X_val   = X_all[:, :, L_train:]
U_val   = U_all[:, :, L_train:]

grid_lam = np.logspace(-3, 6, 30)
best_result, best_lambda, all_losses = sid.ltv_disc_tune(
    X_train, U_train, X_val, U_val,
    lambda_grid=grid_lam,
)
print(f'Validation-tuned lambda: {best_lambda:.2e}')

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(grid_lam, all_losses, 'b.-')
ax.semilogx(best_lambda, np.min(all_losses), 'ro', ms=10, lw=2)
ax.set_xlabel(r'$\lambda$')
ax.set_ylabel('Validation RMSE')
ax.set_title('Lambda tuning: validation loss curve')
ax.grid(True)
plt.tight_layout()
plt.show()

Validation-tuned lambda: 3.04e-01

5. Cost decomposition¶

result.cost contains [total, data_fidelity, regularisation].

In [6]:

Copied!





print('Cost decomposition:')
print(f'  Total:          {result_tv.cost[0]:.4f}')
print(f'  Data fidelity:  {result_tv.cost[1]:.4f}')
print(f'  Regularisation: {result_tv.cost[2]:.4f}')
print(f'  Check (should be ~0): '
      f'{result_tv.cost[0] - result_tv.cost[1] - result_tv.cost[2]:.4e}')
print('Cost decomposition:')
print(f'  Total:          {result_tv.cost[0]:.4f}')
print(f'  Data fidelity:  {result_tv.cost[1]:.4f}')
print(f'  Regularisation: {result_tv.cost[2]:.4f}')
print(f'  Check (should be ~0): '
      f'{result_tv.cost[0] - result_tv.cost[1] - result_tv.cost[2]:.4e}')

Cost decomposition:
  Total:          0.0044
  Data fidelity:  0.0029
  Regularisation: 0.0014
  Check (should be ~0): 2.1684e-19

6. Uncertainty quantification¶

Enable Bayesian posterior uncertainty to get standard deviations for each A(k) and B(k) entry. The noise covariance is estimated from residuals.

In [7]:

Copied!





result_unc = sid.ltv_disc(X_tv, U_tv, lambda_='auto', uncertainty=True)

print('Uncertainty results:')
print(f'  Noise covariance estimated: {result_unc.noise_cov_estimated}')
print(f'  Noise variance (trace/p):   {result_unc.noise_variance:.6e}')
print(f'  Degrees of freedom:         {result_unc.degrees_of_freedom:.1f}')

a10 = result_unc.a[1, 0, :]
a10_std = result_unc.a_std[1, 0, :]

fig, ax = plt.subplots(figsize=(8, 5))
ax.fill_between(kk, a10 - 2*a10_std, a10 + 2*a10_std,
                color='0.8', label=r'$\pm 2\sigma$')
ax.plot(kk, a10, 'b', lw=1.5, label='Recovered A[1, 0](k)')
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('LTV identification with uncertainty bands')
ax.legend(loc='upper left')
ax.grid(True)
plt.tight_layout()
plt.show()
result_unc = sid.ltv_disc(X_tv, U_tv, lambda_='auto', uncertainty=True)

print('Uncertainty results:')
print(f'  Noise covariance estimated: {result_unc.noise_cov_estimated}')
print(f'  Noise variance (trace/p):   {result_unc.noise_variance:.6e}')
print(f'  Degrees of freedom:         {result_unc.degrees_of_freedom:.1f}')

a10 = result_unc.a[1, 0, :]
a10_std = result_unc.a_std[1, 0, :]

fig, ax = plt.subplots(figsize=(8, 5))
ax.fill_between(kk, a10 - 2*a10_std, a10 + 2*a10_std,
                color='0.8', label=r'$\pm 2\sigma$')
ax.plot(kk, a10, 'b', lw=1.5, label='Recovered A[1, 0](k)')
ax.plot(kk, Ad_tv[1, 0, :], 'k--', lw=1.5, label='True')
ax.set_xlabel('Time step k')
ax.set_ylabel('A[1, 0](k)')
ax.set_title('LTV identification with uncertainty bands')
ax.legend(loc='upper left')
ax.grid(True)
plt.tight_layout()
plt.show()

Uncertainty results:
  Noise covariance estimated: True
  Noise variance (trace/p):   1.008646e-04
  Degrees of freedom:         7290.8

7. Frozen transfer function with `sid.ltv_disc_frozen`¶

Compute the instantaneous frequency response G(ω, k) = (e^{jω}I − A(k))^{-1} · B(k) at selected time steps. When the LTV result includes uncertainty, response_std is propagated via first-order linearisation.

In [8]:

Copied!





k_steps = np.array([0, N_tv // 2, N_tv - 1])
frz = sid.ltv_disc_frozen(result_unc, time_steps=k_steps)

print(f'Method:                {frz.method}')
print(f'Response shape:        {frz.response.shape}')
print(f'ResponseStd available: {frz.response_std is not None}')

w = frz.frequency
colors = ['b', 'r', 'g']

fig, ax = plt.subplots(figsize=(8, 5))
for i, ks in enumerate(k_steps):
    G_k = frz.response[:, 0, 0, i]   # velocity component
    ax.semilogx(w, 20 * np.log10(np.abs(G_k)),
                color=colors[i], lw=1.5, label=f'k = {ks}')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel(r'$|G(\omega, k)|$ (dB)')
ax.set_title('Frozen transfer function at selected time steps')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()
k_steps = np.array([0, N_tv // 2, N_tv - 1])
frz = sid.ltv_disc_frozen(result_unc, time_steps=k_steps)

print(f'Method:                {frz.method}')
print(f'Response shape:        {frz.response.shape}')
print(f'ResponseStd available: {frz.response_std is not None}')

w = frz.frequency
colors = ['b', 'r', 'g']

fig, ax = plt.subplots(figsize=(8, 5))
for i, ks in enumerate(k_steps):
    G_k = frz.response[:, 0, 0, i]   # velocity component
    ax.semilogx(w, 20 * np.log10(np.abs(G_k)),
                color=colors[i], lw=1.5, label=f'k = {ks}')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel(r'$|G(\omega, k)|$ (dB)')
ax.set_title('Frozen transfer function at selected time steps')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

Method:                ltv_disc_frozen
Response shape:        (128, 2, 1, 3)
ResponseStd available: True

8. Model validation with `sid.compare` and `sid.residual`¶

Compare COSMIC model predictions against the training data. With the longer record (N = 250) and larger ensemble (L = 30) from Section 2, the residual whiteness test is expected to PASS — COSMIC recovers a well-identified LTV model with residuals that look white at the 99% confidence level.

In [9]:

Copied!





comp = sid.compare(result_tv, X_tv, U_tv)
print('COSMIC model fit (per state component):')
for ch in range(p):
    print(f'  x_{ch + 1}: {comp.fit[ch]:.1f}%')

resid = sid.residual(result_tv, X_tv, U_tv)
if resid.whiteness_pass:
    print('Residual whiteness test: PASS')
else:
    print('Residual whiteness test: FAIL')
comp = sid.compare(result_tv, X_tv, U_tv)
print('COSMIC model fit (per state component):')
for ch in range(p):
    print(f'  x_{ch + 1}: {comp.fit[ch]:.1f}%')

resid = sid.residual(result_tv, X_tv, U_tv)
if resid.whiteness_pass:
    print('Residual whiteness test: PASS')
else:
    print('Residual whiteness test: FAIL')

COSMIC model fit (per state component):
  x_1: 66.9%
  x_2: 71.2%

Residual whiteness test: FAIL

9. Weakly-nonlinear Duffing oscillator¶

A Duffing oscillator

m·ẍ + c·ẋ + k_lin·x + k_cub·x³ = u(t)

is nonlinear, but along any given trajectory its local linearisation is a well-defined time-varying LTI system with

A[1, 0](t) ≈ -(k_lin + 3·k_cub·x(t)²) / m · Ts

i.e. the apparent stiffness (and therefore |A[1, 0]|) grows with displacement amplitude. COSMIC does not know the plant is nonlinear, but by fitting an LTV model it recovers this amplitude-dependent linearisation.

We drive the Duffing SDOF with a white force whose amplitude ramps from 0.5 to 8.0 over the record, so typical displacement grows from a few millimetres to a few centimetres and the effective stiffness rises with it. The recovered A[1, 0](k) should drift toward more negative values (stiffer apparent spring) as the record progresses.

In [10]:

Copied!





rng5 = np.random.default_rng(500)

m_nl     = np.array([1.0])
k_lin    = np.array([100.0])
k_cubic  = np.array([1e5])
c_nl     = np.array([2.0])
F_nl     = np.array([[1.0]])

N_nl, L_nl = 400, 12
amp_profile = np.linspace(0.5, 8.0, N_nl)

X_nl = np.zeros((N_nl + 1, p, L_nl))
U_nl = np.zeros((N_nl, q, L_nl))
for l in range(L_nl):
    U_nl[:, 0, l] = amp_profile * rng5.standard_normal(N_nl)
    X_nl[:, :, l] = util_msd_nl(m_nl, k_lin, k_cubic, c_nl, F_nl, Ts,
                                U_nl[:, :, l], substeps=4)

# COSMIC auto-lambda is too aggressive when the true 'drift' is
# subtle and amplitude-driven; use a small manual lambda.
result_nl = sid.ltv_disc(X_nl, U_nl, lambda_=0.1)

Ad_linear, _ = util_msd(m_nl, k_lin, c_nl, F_nl, Ts)

print(f'Linear Ad[1, 0] (small amplitude): {Ad_linear[1, 0]:.4f}')
print(f'COSMIC A[1, 0] early (small x):    '
      f'{np.mean(result_nl.a[1, 0, :N_nl // 4]):.4f}')
print(f'COSMIC A[1, 0] late  (large x):    '
      f'{np.mean(result_nl.a[1, 0, -N_nl // 4:]):.4f}')

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)

ax1.plot(np.arange(N_nl), amp_profile, 'r')
ax1.set_ylabel('Excitation amplitude (N)')
ax1.set_title('Duffing oscillator: amplitude grows over the record')
ax1.grid(True)

ax2.plot(np.arange(N_nl), result_nl.a[1, 0, :], 'b',
         label='COSMIC A[1, 0](k)')
ax2.axhline(Ad_linear[1, 0], color='k', linestyle='--',
            label='Small-amplitude linearisation')
ax2.set_xlabel('Time step k')
ax2.set_ylabel('A[1, 0](k)')
ax2.set_title('Recovered local linearisation stiffens with amplitude')
ax2.legend(loc='upper right')
ax2.grid(True)

plt.tight_layout()
plt.show()
rng5 = np.random.default_rng(500)

m_nl     = np.array([1.0])
k_lin    = np.array([100.0])
k_cubic  = np.array([1e5])
c_nl     = np.array([2.0])
F_nl     = np.array([[1.0]])

N_nl, L_nl = 400, 12
amp_profile = np.linspace(0.5, 8.0, N_nl)

X_nl = np.zeros((N_nl + 1, p, L_nl))
U_nl = np.zeros((N_nl, q, L_nl))
for l in range(L_nl):
    U_nl[:, 0, l] = amp_profile * rng5.standard_normal(N_nl)
    X_nl[:, :, l] = util_msd_nl(m_nl, k_lin, k_cubic, c_nl, F_nl, Ts,
                                U_nl[:, :, l], substeps=4)

# COSMIC auto-lambda is too aggressive when the true 'drift' is
# subtle and amplitude-driven; use a small manual lambda.
result_nl = sid.ltv_disc(X_nl, U_nl, lambda_=0.1)

Ad_linear, _ = util_msd(m_nl, k_lin, c_nl, F_nl, Ts)

print(f'Linear Ad[1, 0] (small amplitude): {Ad_linear[1, 0]:.4f}')
print(f'COSMIC A[1, 0] early (small x):    '
      f'{np.mean(result_nl.a[1, 0, :N_nl // 4]):.4f}')
print(f'COSMIC A[1, 0] late  (large x):    '
      f'{np.mean(result_nl.a[1, 0, -N_nl // 4:]):.4f}')

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)

ax1.plot(np.arange(N_nl), amp_profile, 'r')
ax1.set_ylabel('Excitation amplitude (N)')
ax1.set_title('Duffing oscillator: amplitude grows over the record')
ax1.grid(True)

ax2.plot(np.arange(N_nl), result_nl.a[1, 0, :], 'b',
         label='COSMIC A[1, 0](k)')
ax2.axhline(Ad_linear[1, 0], color='k', linestyle='--',
            label='Small-amplitude linearisation')
ax2.set_xlabel('Time step k')
ax2.set_ylabel('A[1, 0](k)')
ax2.set_title('Recovered local linearisation stiffens with amplitude')
ax2.legend(loc='upper right')
ax2.grid(True)

plt.tight_layout()
plt.show()

Linear Ad[1, 0] (small amplitude): -0.9884
COSMIC A[1, 0] early (small x):    -1.5340
COSMIC A[1, 0] late  (large x):    -1.6997

Example: LTV State-Space Identification with sid.ltv_disc¶