Frequency-Response Analysis of MIMO System

Open Live Script

This example shows how to estimate the multi-input/multi-output (MIMO) transfer function for a two-body oscillator. This example starts by computing the system inputs and outputs in the time domain, emulating measurements. After applying the definition of the Z-transform to calculate the frequency-domain transfer function of the system, the example uses the functions tfestimate, modalfrf, and modalfit to estimate the transfer function from the simulated measurement data. The comparison between the outputs of these functions shows how you can retrieve the transfer function from system input/output data.

Two-Body Oscillating System

An ideal one-dimensional discrete-time oscillating system consists of two masses, $m_{1}$ and $m_{2}$ , confined between two walls. The units are such that $m_{1} = 1$ kg and $m_{2} = μ$ kg. Each mass is attached to the nearest wall by a spring with elastic constant $k$ (in N/m). An identical spring connects the two masses. Three dampers impede the motion of the masses by exerting on them forces proportional to speed, with damping constant $b$ (in kg/s). Two sensors sample $r_{1}$ and $r_{2}$ (in $m$ ), the displacements of the masses, at a rate $F_{s} = 40$ Hz.

Generate 24,000 time samples, equivalent to 10 minutes. Define the sampling interval $Δ t = 1 / F_{s}$ .

Fs = 40;
dt = 1/Fs;
N = 24000;
t = dt*(0:N-1);

The system can be described by the state-space model

$\begin{array}{c} x (k + 1) = A x (k) + B u (k), \\ y (k) = C x (k) + D u (k), \end{array}$

where $x = {[\begin{array}{c} r_{1} & v_{1} & r_{2} & v_{2} \end{array}]}^{T}$ is the state vector, $r_{i}$ and $v_{i}$ are respectively the location and the velocity of the $i$ th mass, $u = {[\begin{array}{c} u_{1} & u_{2} \end{array}]}^{T}$ is the vector of input driving forces, and $y = {[\begin{array}{c} r_{1} & r_{2} \end{array}]}^{T}$ is the output vector. The state-space matrices are

$A = \exp (A_{c} Δ t), B = A_{c}^{- 1} (A - I) B_{c}, C = [\begin{array}{c} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array}], D = [\begin{array}{c} 0 & 0 \\ 0 & 0 \end{array}],$

$I$ is the $4 \times 4$ identity matrix, and the continuous-time state-space matrices are

$A_{c} = [\begin{array}{c} 0 & 1 & 0 & 0 \\ - 2 k & - 2 b & k & b \\ 0 & 0 & 0 & 1 \\ k / μ & b / μ & - 2 k / μ & - 2 b / μ \end{array}], B_{c} = [\begin{array}{c} 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 / μ \end{array}] .$

Set $k = 400$ N/m, $b = 0.1$ kg/s, and $μ = 0.1$ .

k = 400;
b = 0.1;
mu = 0.1;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/mu b/mu -2*k/mu -2*b/mu];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/mu];
B = Ac\(A-eye(4))*Bc;
C = [1 0 0 0;0 0 1 0];
D = zeros(2);

A random input drives the masses for the first 390 seconds and then the masses are left to rest. Set the initial conditions of location and velocity to zero. Use the state-space model to compute the time evolution of the system. Plot the displacements of the masses as a function of time.

rng("default")
u = randn(2,N);
u(:,t>390) = 0;

y = [0;0];
x = [0;0;0;0];
for iter = 1:N
    y(:,iter) = C*x + D*u(:,iter);
    x = A*x + B*u(:,iter);
end

Set the system time, input, and response data to vertical vector format.

t = t';
u = u';
y = y';

Plot the system response over time for the two masses.

plot3([1 2].*ones(size(t)),t,y)
set(gca,Ydir="reverse")
xticks([1 2])
xticklabels("Mass m_"+[1 2])
ylabel("Time (s)")
zlabel("Displacement (m)")
xlim([0.5 2.5])
grid

Figure contains an axes object. The axes object with ylabel Time (s) contains 2 objects of type line.

Save the system data. The examples Frequency-Response Function of Two-Body Oscillator, Modal Analysis of Two-Body Oscillator, and Transfer Function Estimation of MIMO System use this data file.

save MIMOdata

Theoretical Frequency Response

The frequency-response function of a discrete-time system is the Fourier transform of its time-domain transfer function or, equivalently, the Z-transform of the transfer function evaluated at the unit circle.

Compute the theoretical frequency-response functions. Use 2048 sampling points to calculate the discrete Fourier transform.

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

nfs = 2048;
fz = 0:Fs/nfs:Fs/2;
ztf(1,:,1) = freqz(b1(1,:),a1,fz,Fs);
ztf(1,:,2) = freqz(b1(2,:),a1,fz,Fs);
ztf(2,:,1) = freqz(b2(1,:),a2,fz,Fs);
ztf(2,:,2) = freqz(b2(2,:),a2,fz,Fs);

Transfer Function Estimate using `tfestimate`

Use the system data and the function tfestimate without output arguments to plot the estimate of the MIMO transfer functions. Select the "mimo" option to produce all four transfer functions. Use a periodic 5000-sample Hann window to divide the signals into segments. Specify 2500 samples of overlap between adjoining segments. Store the transfer function of the system as a function of frequency.

wind = hann(5000,"periodic");
nov = 2500;

[tXY,ft] = tfestimate(u,y,wind,nov,nfs,Fs,"mimo");

Verify that the estimate computed by tfestimate coincides with the definition.

plotComparison(fz,ztf,ft,tXY,"tfestimate")

Transfer Function Estimate Using Modal Analysis

Estimate the modal frequency response function of the system using the function modalfrf. Specify the sensor data type as measured displacements.

[frf,f] = modalfrf(u,y,Fs,wind,nov,Sensor="dis");

Plot the estimates and overlay the theoretical predictions.

plotComparison(fz,ztf,f,frf,"modalfrf")

Use modalsd with no output arguments to generate a stabilization diagram for modal analysis and estimate the minimum number of modes. Use a maximum of 10 modes and pick the least-squares rational function estimation method for the calculation.

figure
modalsd(frf,f,Fs,MaxModes=10,FitMethod="lsrf")

Figure contains an axes object. The axes object with title Stabilization Diagram, xlabel Frequency (Hz), ylabel Model Order contains 4 objects of type line. One or more of the lines displays its values using only markers These objects represent Stable in frequency, Stable in frequency and damping, Not stable in frequency, Averaged response function.

The stabilization diagram shows two "+" marks for 2, 4, 5 and higher model order, indicating a stable modal fit for both frequency and damping using at least two modes.

Estimate the natural frequencies, damping ratios, and mode shapes of the system. Use the function modalfit with two modes and pick the least-squares rational function estimation method for the calculation. Obtain the reconstructed transfer functions of the system in the frequency domain.

nModes = 2;
[fn,dr,ms,ofrf] = modalfit(frf,f,Fs,nModes,FitMethod="lsrf");

Compare the natural frequencies to the theoretical predictions for an undamped system.

fn_theo = sqrt(eig([2*k -k;-k/mu 2*k/mu]))/(2*pi);
disp(table(fn,fn_theo,dr,RowNames= "Mode "+(1:nModes)', ...
    VariableNames=["fn" "fn_theo" "dr"]))

                fn      fn_theo       dr    
              ______    _______    _________

    Mode 1    3.8476     3.847     0.0035235
    Mode 2    14.426    14.426      0.011308

The proximity in value of the damping ratio (dr) to zero indicates an underdamped system. The natural frequencies of the system (fn) approach the equivalent undamped natural frequency (fn_theo).

Compare the transfer-function definition with the recovered transfer function from the output of modalfit.

plotComparison(fz,ztf,f,ofrf,"modalfit")

The natural frequencies agree between the recovered and theoretical frequency-response functions despite the differences outside the natural frequencies, correlating with the two modes calculated with modalfit.

Comparison

Compare all transfer function estimation methods discussed in this example.

plotComparison(fz,ztf,ft,tXY,"tfestimate",...
    f,frf,"modalfrf",f,ofrf,"modalfit")

This example uses the functions tfestimate, modalfrf, modalsd, and modalfit to perform frequency-response analysis based on the space-state model of a mass-spring-damper oscillator MIMO system. The functions tfestimate and modalfrf estimate and plot the frequency response of the system, given the information about the system input and output in the time domain. The function modalsd helps to identify the number of modes to use for modal fit. The modal parameters calculated using the function modalfit also help retrieve the frequency response of the system.

Helper Function

The function plotComparison plots and formats a comparison between theoretical and estimated transfer functions.

function plotComparison(varargin)
figure
tiledlayout flow
for jk = 1:2
    for kj = 1:2
        nexttile
        plot(varargin{1}, ...
            mag2db(abs(varargin{2}(jk,:,kj))),"--",LineWidth=2)
        hold on
        for it = 1:fix(nargin/3)
            plot(varargin{3*it},mag2db(abs(varargin{3*it+1}(:,jk,kj))))
        end
        hold off
        grid on
        axis tight
        title("Input "+jk+", Output "+kj)
        xlabel("Frequency (Hz)")
        ylabel("Magnitude (dB)")
    end
    legend("tf definition",varargin{5:3:end},Location="EastOutside")
end
end