Example: Output-COSMIC — LTV Identification from Partial Observations¶
This example demonstrates the Output-COSMIC workflow
(sid.ltv_disc_io) for identifying time-varying systems when
only partial state measurements are available:
$$x(k+1) = A(k)\,x(k) + B(k)\,u(k)$$ $$y(k) = H\,x(k)$$
The workflow is:
- Estimate the frequency response from input-output data.
- Estimate the model order (state dimension) via Hankel SVD.
- Specify the observation matrix
H. - Run
sid.ltv_disc_ioto identify(A, B, x)jointly. - Validate via observation reconstruction (
H · x̂ ≈ y).
Plant. A 2-mass spring-damper chain
wall--k₁--m₁--k₂--m₂ — a natural partial-observation
scenario. The state is [x₁, x₂, v₁, v₂] (two positions and
two velocities), so n = 4. We have position sensors on both
masses but no velocity sensors, so the measurement is
y = [x₁; x₂] (py = 2) and the two velocities are hidden.
Output-COSMIC must infer the hidden velocity trajectories along
with the dynamics.
Gauge note. sid.ltv_disc_io recovers (A, B, x) up to an
unobservable similarity transform — the latent state basis can
be any nonsingular map of the true basis. Direct element-wise
comparison of the recovered A/B against the simulation's
Ad/Bd is therefore not meaningful; what is meaningful is
whether the recovered model predicts the measured outputs. We
check that via the observation reconstruction error at the end.
Translated from exampleOutputCOSMIC.m.
import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd
import sid
%matplotlib inline
System setup¶
The 2-mass SMD gives a 4th-order plant. We pick force actuation at mass 1 and position measurements at both masses.
rng = np.random.default_rng(42)
# ---- Physical plant ----
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]]) # force at mass 1
Ts = 0.01
Ad, Bd = util_msd(m, k, c, F, Ts) # (4, 4), (4, 1)
# ---- Observation matrix: positions only (velocities hidden) ----
n, q, py = 4, 1, 2
H_obs = np.array([[1, 0, 0, 0], # y₁ = x₁
[0, 1, 0, 0]]) # y₂ = x₂
print('Discrete dynamics matrix Ad:')
print(Ad)
print('\nObservation matrix H (measure positions, velocities hidden):')
print(H_obs)
Discrete dynamics matrix Ad: [[ 9.91160615e-01 3.91257699e-03 9.77404408e-03 1.10932852e-04] [ 3.90601777e-03 9.96057895e-01 1.10932852e-04 9.88825654e-03] [-1.75045331e+00 7.73048899e-01 9.52286304e-01 2.32387994e-02] [ 7.71092610e-01 -7.82185895e-01 2.32387994e-02 9.76503248e-01]] Observation matrix H (measure positions, velocities hidden): [[1 0 0 0] [0 1 0 0]]
Simulate trajectories¶
L = 10 trajectories of N = 80 steps each, driven by white
force input with moderate process and measurement noise. The
input is scaled up so that displacements reach a few
centimetres — large enough for the model-order and identification
steps to work cleanly.
N = 80
L = 10
sigma_proc = 1e-3
sigma_meas = 1e-4
X = np.zeros((N + 1, n, L))
U = 5.0 * rng.standard_normal((N, q, L)) # scaled excitation
Y = np.zeros((N + 1, py, L))
for l in range(L):
X[0, :, l] = 0.05 * rng.standard_normal(n)
Y[0, :, l] = H_obs @ X[0, :, l] + sigma_meas * rng.standard_normal(py)
for step in range(N):
X[step + 1, :, l] = (Ad @ X[step, :, l]
+ Bd.ravel() * U[step, 0, l]
+ sigma_proc * rng.standard_normal(n))
Y[step + 1, :, l] = H_obs @ X[step + 1, :, l] + sigma_meas * rng.standard_normal(py)
print(f'Max |Y|: {np.max(np.abs(Y)):.3e} m (measured positions)')
print(f'Max |hidden v|: {np.max(np.abs(X[:, 2:, :])):.3e} m/s')
Max |Y|: 1.061e-01 m (measured positions) Max |hidden v|: 1.053e+00 m/s
Step 1: estimate frequency response¶
Use the first trajectory. Trim Y to match U's length
(sid.freq_bt requires equal-length signals).
G = sid.freq_bt(Y[:N, :, 0], U[:, :, 0], window_size=20, sample_time=Ts)
print(f'G.response shape: {G.response.shape}')
G.response shape: (128, 2, 1)
Step 2: model-order determination¶
sid.model_order estimates the state dimension from Hankel
singular values. For short records with sharp resonances the
estimate can overshoot the true order; we clip to the known
n = 4 to keep the downstream dimensions consistent.
n_est, sv_info = sid.model_order(G)
print(f'Hankel-SVD estimate: n = {n_est} (true = {n})')
print('(For short lightly-damped records the Hankel estimate can overshoot;')
print(" we pin n = 4 for the identification step below.)")
Hankel-SVD estimate: n = 19 (true = 4) (For short lightly-damped records the Hankel estimate can overshoot; we pin n = 4 for the identification step below.)
Step 3: construct observation matrix¶
In practice H encodes the known sensor layout. Here we pass
the true position-measurement matrix since that's what we built
into the simulation.
H_use = H_obs.copy()
print(f'Observation matrix H ({H_use.shape[0]} x {H_use.shape[1]}):')
print(H_use)
Observation matrix H (2 x 4): [[1 0 0 0] [0 1 0 0]]
Step 4: identify the LTV model via sid.ltv_disc_io¶
Runs Output-COSMIC: jointly infers the latent state trajectory
and the discrete (A(k), B(k)) consistent with both the
dynamics and the partial observations.
print('Running Output-COSMIC identification...')
result = sid.ltv_disc_io(Y, U, H_use, lambda_=1e5)
print(f'Converged in {result.iterations} iterations.')
print(f'Final cost: {result.cost[-1]:.4f}')
Running Output-COSMIC identification...
Converged in 22 iterations. Final cost: 0.0052
Convergence history¶
fig, ax = plt.subplots(figsize=(8, 4))
ax.semilogy(np.arange(1, len(result.cost) + 1), result.cost, 'b-o', markersize=4)
ax.set_xlabel('Iteration')
ax.set_ylabel('Cost J')
ax.set_title('Output-COSMIC: convergence')
ax.grid(True)
plt.tight_layout()
plt.show()
State recovery: observed channels vs hidden channels¶
Position channels (measured) should match the true positions
closely. Velocity channels (hidden — never measured) are
inferred by Output-COSMIC from the position dynamics; the
inferred curve should track the true velocity up to the
algorithm's gauge transform (here the basis is pinned by H so
the comparison is meaningful).
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
t_axis = np.arange(N + 1)
# Observed channel: position of mass 1
axes[0, 0].plot(t_axis, X[:, 0, 0], 'k-', lw=1.5, label='True x₁')
axes[0, 0].plot(t_axis, result.x[:, 0, 0], 'b--', label='Estimated x₁')
axes[0, 0].plot(t_axis, Y[:, 0, 0], 'r.', ms=4, label='Measurement y₁')
axes[0, 0].legend()
axes[0, 0].set_xlabel('k')
axes[0, 0].set_ylabel('x₁ (m)')
axes[0, 0].set_title('Observed: position of mass 1')
axes[0, 0].grid(True)
# Observed channel: position of mass 2
axes[0, 1].plot(t_axis, X[:, 1, 0], 'k-', lw=1.5, label='True x₂')
axes[0, 1].plot(t_axis, result.x[:, 1, 0], 'b--', label='Estimated x₂')
axes[0, 1].plot(t_axis, Y[:, 1, 0], 'r.', ms=4, label='Measurement y₂')
axes[0, 1].legend()
axes[0, 1].set_xlabel('k')
axes[0, 1].set_ylabel('x₂ (m)')
axes[0, 1].set_title('Observed: position of mass 2')
axes[0, 1].grid(True)
# Hidden channel: velocity of mass 1
axes[1, 0].plot(t_axis, X[:, 2, 0], 'k-', lw=1.5, label='True v₁')
axes[1, 0].plot(t_axis, result.x[:, 2, 0], 'b--', label='Estimated v₁')
axes[1, 0].legend()
axes[1, 0].set_xlabel('k')
axes[1, 0].set_ylabel('v₁ (m/s)')
axes[1, 0].set_title('Hidden: velocity of mass 1')
axes[1, 0].grid(True)
# Hidden channel: velocity of mass 2
axes[1, 1].plot(t_axis, X[:, 3, 0], 'k-', lw=1.5, label='True v₂')
axes[1, 1].plot(t_axis, result.x[:, 3, 0], 'b--', label='Estimated v₂')
axes[1, 1].legend()
axes[1, 1].set_xlabel('k')
axes[1, 1].set_ylabel('v₂ (m/s)')
axes[1, 1].set_title('Hidden: velocity of mass 2')
axes[1, 1].grid(True)
plt.tight_layout()
plt.show()
Validation: observation reconstruction error¶
The gauge-invariant test. Compute H · x̂ from the recovered
state and compare against the measured y. This metric is
well-defined regardless of any similarity transform that
Output-COSMIC's basis may have applied internally.
Y_recon = np.zeros((N + 1, py, L))
for l in range(L):
Y_recon[:, :, l] = result.x[:, :, l] @ H_use.T
obs_err = np.linalg.norm(Y_recon.ravel() - Y.ravel()) / np.linalg.norm(Y.ravel())
print(f'Observation reconstruction error: {obs_err:.4f}')
print(f'(Relative Frobenius over all {N+1} time steps x {py} channels x {L} trajectories.)')
Observation reconstruction error: 0.0174 (Relative Frobenius over all 81 time steps x 2 channels x 10 trajectories.)
Frozen-time inspection of the recovered A and B¶
Print the recovered A(k) and B(k) at the midpoint of the
record. Remember these are in a latent basis — direct
comparison to Ad and Bd is not meaningful, but you should
still see sensible magnitudes (entries O(1), not blowing up).
mid_k = N // 2
print(f'A at midpoint (k = {mid_k}):')
print(result.a[:, :, mid_k])
print(f'\nB at midpoint (k = {mid_k}):')
print(result.b[:, :, mid_k])
A at midpoint (k = 40): [[ 1.73693923e+00 2.30147418e-01 2.59112899e-03 -2.28344182e-02] [ 1.78791708e-01 1.01913329e+00 2.07521088e-03 -7.17683736e-03] [ 5.05813600e+01 1.39452866e+01 3.49634968e-01 -1.59624103e-01] [ 3.63158920e+01 1.04766171e+01 4.11408202e-01 -5.60485487e-01]] B at midpoint (k = 40): [[ 1.85704492e-05] [-2.36546660e-05] [ 5.76549931e-03] [-4.27424177e-03]]