基于特征系统实现算法 (ERA) 的系统辨识

此示例使用:

本示例考虑了一个由一系列锤击激励的 3 自由度 (DOF) 系统。产生的位移由传感器记录下来。该系统采用比例阻尼，因此阻尼矩阵是质量矩阵和刚度矩阵的线性组合。本示例中，模态分析部分需要使用 Signal Processing Toolbox™。

预处理数据

导入第一组测量数据。这些数据包括激励信号、响应信号、时间信号以及真实频率响应函数。响应信号（记为 Y1）表示第一个质量的位移。激励信号由十次连续的锤击组成，而响应信号包含相应的位移。每个撞击信号的持续时间为 2.53 秒。激励信号和响应信号中混入了加性噪声。

load modaldata XhammerMISO1 YhammerMISO1 fs; 
rng('default');

% Add noise to excitation and response signals
XhammerMISO1 = XhammerMISO1 + randn(size(YhammerMISO1))/1250;
YhammerMISO1 = YhammerMISO1 + randn(size(YhammerMISO1))/1e11;

% Define outputs from function
t = (0:size(XhammerMISO1,1)-1)/fs';
X1 = 1e2*XhammerMISO1;
Y1 = 1e2*YhammerMISO1;
X0 = X1(:,1);
Y0 = Y1(:,1);

可视化测量的第一个激励通道和响应通道。

subplot(2,1,1)
plot(t,X0)
xlabel('Time (s)')
ylabel('Force (N)')
grid on;
title('Excitation and Response for a 3DOF System')
subplot(2,1,2)
plot(t,Y0)
xlabel('Time (s)')
ylabel('Displacement (m)')
grid on;

Figure contains 2 axes objects. Axes object 1 with title Excitation and Response for a 3DOF System, xlabel Time (s), ylabel Force (N) contains an object of type line. Axes object 2 with xlabel Time (s), ylabel Displacement (m) contains an object of type line.

使用 ERA 辨识系统

为估计数据和验证数据创建 iddata 对象。

Ts = 1/fs; % sample time
estimationData = iddata(Y0(1:1000), X0(1:1000), 1/fs);
validationData = iddata(Y0(1001:2000), X0(1001:2000), 1/fs);

对估计数据进行可视化。

figure
plot(estimationData)

输入数据的绘图显示存在输入延迟。从数据中去除延迟。

[~,inputDelay] = max(estimationData.InputData);
estimationData = estimationData(inputDelay:end);

era 函数要求以 timetable 或数值矩阵的形式传递数据。将 iddata 对象转换为 timetable。

L = length(estimationData.InputData);
t = seconds(estimationData.Tstart + (0:L-1)'*Ts);
y1 = estimationData.OutputData;
tt = timetable(t,y1);

使用 era 根据这些估计数据对状态空间模型进行估计。

order = 6;
sys = era(tt, order, 'Feedthrough', true)

sys =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
               x1          x2          x3          x4          x5          x6
   x1      0.8339     -0.5523   -0.001099    0.003053   -0.001497   0.0004633
   x2      0.5523      0.8323    -0.00714   -0.001638   -0.002676   0.0008873
   x3   -0.001099     0.00714      0.2293      0.9709  -0.0009399  -0.0006184
   x4   -0.003053   -0.001638     -0.9709      0.2289    0.006094   -0.002101
   x5   -0.001497    0.002676  -0.0009399   -0.006094     -0.5463     -0.8302
   x6  -0.0004633   0.0008873   0.0006184   -0.002101      0.8302     -0.5465
 
  B = 
               u1
   x1   -0.001403
   x2   0.0007461
   x3   -0.001515
   x4  -0.0001808
   x5   -0.000818
   x6  -0.0002384
 
  C = 
               x1          x2          x3          x4          x5          x6
   y1  -3.508e-07  -1.865e-07  -3.787e-07    4.52e-08  -2.045e-07   5.959e-08
 
  D = 
              u1
   y1  1.852e-10
 
  K = 
       y1
   x1   0
   x2   0
   x3   0
   x4   0
   x5   0
   x6   0
 
Sample time: 0.00025 seconds

Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: yes
   Disturbance component: none
   Number of free coefficients: 49
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                              
Estimated using the Eigensystem Realization Algorithm
 
Model Properties

将估计系统的响应与验证数据进行比较。

compare(validationData,sys);

MATLAB figure

将所估计模型的性能与另一个使用 ssest 估计的状态空间模型进行比较。

sys2 = ssest(estimationData, order, 'Feedthrough', true)

sys2 =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
           x1      x2      x3      x4      x5      x6
   x1   4.985    2323  -66.74   241.3   35.31  -5.076
   x2   -2341  -7.766  -300.1   11.19    11.2   25.28
   x3   49.75     311  -5.616   -5360   -30.9   -15.2
   x4  -270.4  -6.758    5333  -15.06   1.474   64.71
   x5  -52.95  -16.58  -14.17  -34.15   3.021    8548
   x6  -41.44  -15.95   -16.8  -134.7   -8674  -51.74
 
  B = 
              u1
   x1  -0.001651
   x2   -0.02103
   x3    -0.1236
   x4    -0.2095
   x5     -1.006
   x6     0.9286
 
  C = 
               x1          x2          x3          x4          x5          x6
   y1  -3.244e-05   2.517e-05   1.609e-05   1.559e-05   2.731e-06   2.401e-06
 
  D = 
              u1
   y1  5.259e-09
 
  K = 
               y1
   x1  -9.881e+06
   x2   8.603e+06
   x3   4.435e+06
   x4   2.032e+07
   x5   7.865e+07
   x6  -4.392e+07
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: yes
   Disturbance component: estimate
   Number of free coefficients: 55
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                    
Estimated using SSEST on time domain data "estimationData".
Fit to estimation data: 99.6% (prediction focus)           
FPE: 5.214e-17, MSE: 4.898e-17                             
 
Model Properties

compare(validationData,sys,sys2);

MATLAB figure

legend('Validation data','ERA','SSEST');

MATLAB figure

这两个模型的拟合相近。

模态分析

使用来自 Signal Processing Toolbox 包的 modalfit 函数，对使用 era 估计的状态空间模型进行模态分析。modalfit 返回阻尼比 zeta 和阻尼固有频率 fd。

[~,f] = modalfrf(sys);
[fd, zeta] = modalfit(sys,f,3);

查看 zeta。

zeta

zeta = 3×1

    0.0008
    0.0018
    0.0028

查看 fd。

fd

将这些结果与使用 ssest 估计的模型所得的模态进行比较。

[frf,f] = modalfrf(sys2);
[fd, zeta] = modalfit(sys2,f,3)

zeta = 3×1

    0.0008
    0.0018
    0.0029

从 modalfit 方法获得的第二个估计模型的模态与使用 era 方法估计的模型所获得的模态之间非常吻合。

降阶建模比较

比较在系统辨识中使用较低模型阶数时，era 和 ssest 的性能。

order = 4;
sys = era(tt, order, 'Feedthrough', true);
sys2 = ssest(estimationData, order, 'Feedthrough', true);

像之前那样比较这两个模型的性能。

compare(validationData,sys,sys2);
legend('Validation data','ERA','SSEST');

MATLAB figure

对于降阶模型，ERA 模型的拟合效果明显优于 SSEST 模型。该结果证明了 ERA 方法在从带噪声的冲激响应数据中获取状态空间模型方面的能力。

另请参阅

era | iddata | timetable | modalfit (Signal Processing Toolbox) | modalfrf (Signal Processing Toolbox)