💥💥💞💞 欢迎来到本博客 ❤️❤️💥💥
🏆 博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。
⛳️ 座右铭:行百里者,半于九十。
📋📋📋 本文目录如下: 🎁🎁🎁
目录
💥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,这些向量构成了系统的状态空间基。
- 对观测矩阵 Y1Y1 进行奇异值分解,得到:
-
系统矩阵估计:
- 利用 UnUn 和 Y2Y2 计算系统矩阵 AA 和 CC:
A^=Un+Y2VnΣn−1A^=Un+Y2VnΣn−1
C^=Y1VnΣn−1C^=Y1VnΣn−1
- 其中 Un+Un+ 表示 UnUn 的伪逆。
- 利用 UnUn 和 Y2Y2 计算系统矩阵 AA 和 CC:
-
模态参数提取:
- 通过对系统矩阵 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]);
🎉3 参考文献
文章中一些内容引自网络,会注明出处或引用为参考文献,难免有未尽之处,如有不妥,请随时联系删除。
[1]赵欢欢,熊炜,张超,等.基于数字孪生的MMC半桥子模块模型构建及参数辨识[J/OL].现代电力,1-9[2024-11-03].https://doi.org/10.19725/j.cnki.1007-2322.2023.0408.
[2]贾科,刘芸,毕天姝,等.构网型控制的新能源电源故障控制参数辨识方法[J/OL].中国电机工程学报,1-10[2024-11-03].http://kns.cnki.net/kcms/detail/11.2107.tm.20241028.1422.017.html.