【参数辨识】随机子空间辨识 (SSI)（Matlab实现）-优快云博客

本文链接：https://blog.youkuaiyun.com/Matlab88888/article/details/143464141

💥💥💞💞 欢迎来到本博客 ❤️❤️💥💥

🏆 博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️ 座右铭：行百里者，半于九十。

📋📋📋 本文目录如下： 🎁🎁🎁

💥1 概述

随机子空间辨识（Stochastic Subspace Identification， SSI）是一种用于线性动态系统辨识的强大方法。该方法最初由 Peter Van Overschee 和 Bart De Moor 于 1996 年提出，后来在结构健康监测、振动分析、控制系统建模等领域得到了广泛应用。SSI 的主要优点在于其能够从输入-输出数据中高效地估计系统的模态参数（如自然频率、阻尼比和模态形状），并且具有良好的鲁棒性和数值稳定性。

基本原理：

SSI 方法的核心思想是利用系统的输入-输出数据构建一个扩展的观测矩阵，然后通过奇异值分解（SVD）或特征值分解（EVD）技术提取系统的动态特性。具体步骤如下：

数据预处理：
- 收集系统的输入和输出数据。
- 对数据进行预处理，如去除均值、归一化等，以提高后续计算的准确性。
构造观测矩阵：
- 从输入-输出数据中构造两个观测矩阵 Y1Y1 和 Y2Y2，其中 Y1Y1 包含前 N−1N−1 个时间步的数据，Y2Y2 包含后 N−1N−1 个时间步的数据。
- 观测矩阵的形式通常为：
  Y1=[y（1）y（2）⋯y（N−1）y（2）y（3）⋯y（N）⋮⋮⋱⋮y（p）y（p+1）⋯y（p+N−2）]Y1=y（1）y（2）⋮y（p）y（2）y（3）⋮y（p+1）⋯⋯⋱⋯y（N−1）y（N）⋮y（p+N−2）
  
  Y2=[y（2）y（3）⋯y（N）y（3）y（4）⋯y（N+1）⋮⋮⋱⋮y（p+1）y（p+2）⋯y（p+N−1）]Y2=y（2）y（3）⋮y（p+1）y（3）y（4）⋮y（p+2）⋯⋯⋱⋯y（N）y（N+1）⋮y（p+N−1）
奇异值分解（SVD）：
- 对观测矩阵 Y1Y1 进行奇异值分解，得到：
  Y1=UΣVTY1=UΣVT
- 选择前 nn 个最大的奇异值对应的左奇异向量组成矩阵 UnUn，这些向量构成了系统的状态空间基。
系统矩阵估计：
- 利用 UnUn 和 Y2Y2 计算系统矩阵 AA 和 CC：
  A^=Un+Y2VnΣn−1A^=Un+Y2VnΣn−1
  
  C^=Y1VnΣn−1C^=Y1VnΣn−1
- 其中 Un+Un+ 表示 UnUn 的伪逆。
模态参数提取：
- 通过对系统矩阵 A^A^ 进行特征值分解，可以得到系统的模态参数，如自然频率和阻尼比。
- 模态形状可以通过特征向量和矩阵 C^C^ 进行计算。

📚2 运行结果

主函数部分代码：

clc; clear; close all;
%Model Parameters and excitation
%--------------------------------------------------------------------------

M=[1 0; 0 1];
K=[2 -1; -1 1]*5;
C=0.1*M+0.1*K;
f=2*randn(2,10000);
fs=100;

%Apply modal superposition to get response
%--------------------------------------------------------------------------

n=size(f,1);
dt=1/fs; %sampling rate
[Vectors, Values]=eig(K,M);
Freq=sqrt(diag(Values))/(2*pi); % undamped natural frequency
steps=size(f,2);

Mn=diag(Vectors'*M*Vectors); % uncoupled mass
Cn=diag(Vectors'*C*Vectors); % uncoupled damping
Kn=diag(Vectors'*K*Vectors); % uncoupled stifness
wn=sqrt(diag(Values));
zeta=Cn./(sqrt(2.*Mn.*Kn));  % damping ratio
wd=wn.*sqrt(1-zeta.^2);

fn=Vectors'*f; % generalized input force matrix

t=[0:dt:dt*steps-dt];

for i=1:1:n
    
    h(i,:)=(1/(Mn(i)*wd(i))).*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t); %transfer function of displacement
    hd(i,:)=(1/(Mn(i)*wd(i))).*(-zeta(i).*wn(i).*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)+wd(i).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t)); %transfer function of velocity
    hdd(i,:)=(1/(Mn(i)*wd(i))).*((zeta(i).*wn(i))^2.*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)-zeta(i).*wn(i).*wd(i).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t)-wd(i).*((zeta(i).*wn(i)).*exp(-zeta(i)*wn(i)*t).*cos(wd(i)*t))-wd(i)^2.*exp(-zeta(i)*wn(i)*t).*sin(wd(i)*t)); %transfer function of acceleration
    
    qq=conv(fn(i,:),h(i,:))*dt;
    qqd=conv(fn(i,:),hd(i,:))*dt;
    qqdd=conv(fn(i,:),hdd(i,:))*dt;
    
    q(i,:)=qq(1:steps); % modal displacement
    qd(i,:)=qqd(1:steps); % modal velocity
    qdd(i,:)=qqdd(1:steps); % modal acceleration
       
end

x=Vectors*q; %displacement
v=Vectors*qd; %vecloity
a=Vectors*qdd; %vecloity

%Add noise to excitation and response
%--------------------------------------------------------------------------
f2=f+0.05*randn(2,10000);
a2=a+0.05*randn(2,10000);
v2=v+0.05*randn(2,10000);
x2=x+0.05*randn(2,10000);

%Plot displacement of first floor without and with noise
%--------------------------------------------------------------------------
figure;
subplot(3,2,1)
plot(t,f(1,:)); xlabel('Time (sec)');  ylabel('Force1'); title('First Floor');
subplot(3,2,2)
plot(t,f(2,:)); xlabel('Time (sec)');  ylabel('Force2'); title('Second Floor');
subplot(3,2,3)
plot(t,x(1,:)); xlabel('Time (sec)');  ylabel('DSP1');
subplot(3,2,4)
plot(t,x(2,:)); xlabel('Time (sec)');  ylabel('DSP2');
subplot(3,2,5)
plot(t,x2(1,:)); xlabel('Time (sec)');  ylabel('DSP1+Noise');
subplot(3,2,6)
plot(t,x2(2,:)); xlabel('Time (sec)');  ylabel('DSP2+Noise');

%Identify modal parameters using displacement with added uncertainty
%--------------------------------------------------------------------------
output=x2;
ncols=5000;
nrows=600;
cut=4;  %Identify 2 modes
[Result]=SSID(output,fs,ncols,nrows,cut);    %SSI

%Plot real and identified first modes to compare between them
%--------------------------------------------------------------------------
figure;
plot([0 ; -Vectors(:,1)],[0 1 2],'r*-');
hold on
plot([0  ;Result.Parameters.ModeShape(:,1)],[0 1 2],'go-.');
hold on
plot([0 ; -Vectors(:,2)],[0 1 2],'b^--');
hold on
plot([0  ;Result.Parameters.ModeShape(:,2)],[0 1 2],'mv-');
hold off
title('Real and Identified Mode Shapes');
legend('Mode 1 (Real)','Mode 1 (Identified using SSI)','Mode 2 (Real)','Mode 2 (Identified using SSI)');
xlabel('Amplitude');
ylabel('Floor');
grid on;
daspect([1 1 1]);