Example: Multi-Input Multi-Output (MIMO) Frequency Response Estimation¶

This example demonstrates MIMO system identification with sid.freq_bt. Key differences from SISO: response is a 3-D array shaped (nf, ny, nu), noise_spectrum is a Hermitian spectral matrix, and coherence is not available.

Plant. A 2-mass spring-damper chain wall--k₁--m₁--k₂--m₂, naturally a MIMO system: force actuators at each mass form the inputs and position sensors at each mass form the outputs, so we can study full 2×2 transfer matrices.

Parameter	Value
m	[1, 1] kg
k	[100, 80] N/m
c	[2, 2] N·s/m

The two normal modes sit at approximately 6.0 and 15.0 rad/s. Cross-coupling G₁₂ and G₂₁ comes from the physical connection through spring k₂.

Translated from exampleMIMO.m.

In [1]:

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

%matplotlib inline

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

2-output, 1-input system¶

First a simple case: a single force is applied at mass 1, and we measure positions of both masses. This is a 1×2 MIMO response — response.shape == (nf, 2, 1). The static gains should match the first column of inv(K) (the static compliance matrix): [1/k₁, 1/k₁] = [0.01, 0.01] approximately, because the first column of K^{-1} equals [1/k₁, 1/k₁] for this topology.

In [2]:

rng = np.random.default_rng(10)

# ---- 2-mass plant, force only at mass 1 ----
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 only
Ts = 0.01
N  = 4000

Ad, Bd = util_msd(m, k, c, F, Ts)
C_out = np.array([[1.0, 0.0, 0.0, 0.0],
                  [0.0, 1.0, 0.0, 0.0]])   # measure x₁ and x₂

u = rng.standard_normal(N)
x = np.zeros((N + 1, 4))
for step in range(N):
    x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step]
y = x[1:, :2] + 2e-4 * rng.standard_normal((N, 2))   # two position outputs

w_grid = np.linspace(0.005, np.pi, 512)
result = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)

rng = np.random.default_rng(10) # ---- 2-mass plant, force only at mass 1 ---- 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 only Ts = 0.01 N = 4000 Ad, Bd = util_msd(m, k, c, F, Ts) C_out = np.array([[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0]]) # measure x₁ and x₂ u = rng.standard_normal(N) x = np.zeros((N + 1, 4)) for step in range(N): x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step] y = x[1:, :2] + 2e-4 * rng.standard_normal((N, 2)) # two position outputs w_grid = np.linspace(0.005, np.pi, 512) result = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)

Inspect MIMO result dimensions¶

In [3]:

print(f'Response shape:       {result.response.shape}  # (nf, ny, nu)')
print(f'Noise spectrum shape: {result.noise_spectrum.shape}  # (nf, ny, ny)')
print(f'Coherence is None:    {result.coherence is None}')

print(f'Response shape: {result.response.shape} # (nf, ny, nu)') print(f'Noise spectrum shape: {result.noise_spectrum.shape} # (nf, ny, ny)') print(f'Coherence is None: {result.coherence is None}')

Response shape:       (512, 2, 1)  # (nf, ny, nu)
Noise spectrum shape: (512, 2, 2)  # (nf, ny, ny)
Coherence is None:    True

Plot individual output channels¶

sid.bode_plot only shows the first channel, so we plot both manually. The dashed reference is the true discrete transfer function C_y · (e^{jω}·I − Ad)^{-1} · Bd.

In [4]:

w = result.frequency

def true_tf_mimo(w, C_sel, Bd_sel):
    '''Return (nf,) complex transfer from Bd_sel input to C_sel output.'''
    return np.array([
        (C_sel @ np.linalg.inv(np.exp(1j * wi) * np.eye(4) - Ad) @ Bd_sel)[0, 0]
        for wi in w
    ])

G1_true = true_tf_mimo(w, C_out[0:1, :], Bd[:, [0]])
G2_true = true_tf_mimo(w, C_out[1:2, :], Bd[:, [0]])

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

ax1.semilogx(w, 20 * np.log10(np.abs(result.response[:, 0, 0])),
             'b', label='Estimated')
ax1.semilogx(w, 20 * np.log10(np.abs(G1_true)),
             'k--', label='True')
ax1.set_ylabel('Magnitude (dB)')
ax1.set_title(r'$G_{1}$: force at mass 1 $\to$ position $x_1$')
ax1.legend()
ax1.grid(True)

ax2.semilogx(w, 20 * np.log10(np.abs(result.response[:, 1, 0])),
             'r', label='Estimated')
ax2.semilogx(w, 20 * np.log10(np.abs(G2_true)),
             'k--', label='True')
ax2.set_ylabel('Magnitude (dB)')
ax2.set_xlabel('Frequency (rad/sample)')
ax2.set_title(r'$G_{2}$: force at mass 1 $\to$ position $x_2$')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

w = result.frequency def true_tf_mimo(w, C_sel, Bd_sel): '''Return (nf,) complex transfer from Bd_sel input to C_sel output.''' return np.array([ (C_sel @ np.linalg.inv(np.exp(1j * wi) * np.eye(4) - Ad) @ Bd_sel)[0, 0] for wi in w ]) G1_true = true_tf_mimo(w, C_out[0:1, :], Bd[:, [0]]) G2_true = true_tf_mimo(w, C_out[1:2, :], Bd[:, [0]]) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True) ax1.semilogx(w, 20 * np.log10(np.abs(result.response[:, 0, 0])), 'b', label='Estimated') ax1.semilogx(w, 20 * np.log10(np.abs(G1_true)), 'k--', label='True') ax1.set_ylabel('Magnitude (dB)') ax1.set_title(r'$G_{1}$: force at mass 1 $\to$ position $x_1$') ax1.legend() ax1.grid(True) ax2.semilogx(w, 20 * np.log10(np.abs(result.response[:, 1, 0])), 'r', label='Estimated') ax2.semilogx(w, 20 * np.log10(np.abs(G2_true)), 'k--', label='True') ax2.set_ylabel('Magnitude (dB)') ax2.set_xlabel('Frequency (rad/sample)') ax2.set_title(r'$G_{2}$: force at mass 1 $\to$ position $x_2$') ax2.legend() ax2.grid(True) plt.tight_layout() plt.show()

Noise spectral matrix¶

For a 2-output system, noise_spectrum is (nf, 2, 2): a Hermitian positive semi-definite matrix at each frequency. The diagonal elements are the individual output noise spectra; the off-diagonal elements are cross-spectra that encode how unmeasured disturbances correlate between the two outputs.

In [5]:

diag11 = np.real(result.noise_spectrum[:, 0, 0])
diag22 = np.real(result.noise_spectrum[:, 1, 1])

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w, 10 * np.log10(diag11), 'b', label=r'$\Phi_{v,11}$')
ax.semilogx(w, 10 * np.log10(diag22), 'r', label=r'$\Phi_{v,22}$')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Noise spectrum (dB)')
ax.set_title('Diagonal elements of the noise spectral matrix')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

diag11 = np.real(result.noise_spectrum[:, 0, 0]) diag22 = np.real(result.noise_spectrum[:, 1, 1]) fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(w, 10 * np.log10(diag11), 'b', label=r'$\Phi_{v,11}$') ax.semilogx(w, 10 * np.log10(diag22), 'r', label=r'$\Phi_{v,22}$') ax.set_xlabel('Frequency (rad/sample)') ax.set_ylabel('Noise spectrum (dB)') ax.set_title('Diagonal elements of the noise spectral matrix') ax.legend() ax.grid(True) plt.tight_layout() plt.show()

2-output, 2-input system¶

Now a full 2×2 MIMO: independent white forces at both masses, measure both positions. The cross-coupling through spring k₂ gives non-trivial off-diagonal transfers G₁₂ and G₂₁.

util_msd handles multi-input plants directly — pass F = I₂ and the returned Bd has shape (4, 2).

In [6]:

rng2 = np.random.default_rng(20)

F_22 = np.array([[1.0, 0.0],
                 [0.0, 1.0]])   # force actuators at both masses
Ad_22, Bd_22 = util_msd(m, k, c, F_22, Ts)

N2 = 4000
u2in = rng2.standard_normal((N2, 2))
x2 = np.zeros((N2 + 1, 4))
for step in range(N2):
    x2[step + 1] = Ad_22 @ x2[step] + Bd_22 @ u2in[step]
y2out = x2[1:, :2] + 2e-4 * rng2.standard_normal((N2, 2))

result_22 = sid.freq_bt(y2out, u2in, window_size=200,
                        frequencies=w_grid, sample_time=Ts)

print(f'2x2 MIMO response shape: {result_22.response.shape}')

rng2 = np.random.default_rng(20) F_22 = np.array([[1.0, 0.0], [0.0, 1.0]]) # force actuators at both masses Ad_22, Bd_22 = util_msd(m, k, c, F_22, Ts) N2 = 4000 u2in = rng2.standard_normal((N2, 2)) x2 = np.zeros((N2 + 1, 4)) for step in range(N2): x2[step + 1] = Ad_22 @ x2[step] + Bd_22 @ u2in[step] y2out = x2[1:, :2] + 2e-4 * rng2.standard_normal((N2, 2)) result_22 = sid.freq_bt(y2out, u2in, window_size=200, frequencies=w_grid, sample_time=Ts) print(f'2x2 MIMO response shape: {result_22.response.shape}')

2x2 MIMO response shape: (512, 2, 2)

Plot the full 2×2 transfer matrix¶

The diagonal elements G₁₁ and G₂₂ are the direct compliances of each mass (force at mass i, position of mass i). The off-diagonal elements G₁₂ and G₂₁ are the transmission responses through the coupling spring — both exhibit the same two resonances because all modes couple both masses.

In [7]:

titles = [[r'$G_{11}$: $u_1 \to y_1$', r'$G_{12}$: $u_2 \to y_1$'],
          [r'$G_{21}$: $u_1 \to y_2$', r'$G_{22}$: $u_2 \to y_2$']]

fig, axes = plt.subplots(2, 2, figsize=(10, 7), sharex=True)
for iy in range(2):
    for iu in range(2):
        ax = axes[iy, iu]
        ax.semilogx(w, 20 * np.log10(np.abs(result_22.response[:, iy, iu])), 'b')
        ax.set_ylabel('Magnitude (dB)')
        if iy == 1:
            ax.set_xlabel('Frequency (rad/sample)')
        ax.set_title(titles[iy][iu])
        ax.grid(True)

plt.tight_layout()
plt.show()

titles = [[r'$G_{11}$: $u_1 \to y_1$', r'$G_{12}$: $u_2 \to y_1$'], [r'$G_{21}$: $u_1 \to y_2$', r'$G_{22}$: $u_2 \to y_2$']] fig, axes = plt.subplots(2, 2, figsize=(10, 7), sharex=True) for iy in range(2): for iu in range(2): ax = axes[iy, iu] ax.semilogx(w, 20 * np.log10(np.abs(result_22.response[:, iy, iu])), 'b') ax.set_ylabel('Magnitude (dB)') if iy == 1: ax.set_xlabel('Frequency (rad/sample)') ax.set_title(titles[iy][iu]) ax.grid(True) plt.tight_layout() plt.show()

MIMO uncertainty¶

In v1.0, response_std is NaN for MIMO (no asymptotic formula implemented).

In [8]:

print(f'MIMO response_std contains NaN: {np.all(np.isnan(result_22.response_std))}')

print(f'MIMO response_std contains NaN: {np.all(np.isnan(result_22.response_std))}')

MIMO response_std contains NaN: False