使用主成分分析进行模态分解(Matlab代码实现)

 👨‍🎓个人主页:研学社的博客 

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

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

⛳️座右铭:行百里者,半于九十。

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

目录

💥1 概述

📚2 运行结果

🌈3 Matlab代码实现

🎉4 参考文献


💥1 概述

本文使用主成分分析 (PCA) 对 2DOF 系统进行振型识别,该系统受到高斯白噪声激励,并增加了响应的不确定性(也是高斯白噪声)。但是,由于
协方差矩阵的对称性,PCA的特征向量是正颌的。
-模态形状仅在矩阵 inv(M)*K 对称时才是正交的。
-PCA 只有在实振型是正颌时才识别它们,这意味着 inv(M)*K 是对称的。

通过改变质量矩阵M=[2 0; 0 1];而不是恒等式inv(M)*K将不是对称的,即使刚度矩阵是对称的,PCA将无法识别实际的振型。

📚2 运行结果

 

🌈3 Matlab代码实现

部分代码:

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
%--------------------------------------------------------------------------
[V]=pca(x2');   %PCA eigenvectors

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

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
%--------------------------------------------------------------------------
[V]=pca(x2');   %PCA eigenvectors

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

🎉4 参考文献

部分理论来源于网络,如有侵权请联系删除。

[1] Al Rumaithi, Ayad, "Characterization of Dynamic Structures Using Parametric and Non-parametric System Identification Methods" (2014). Electronic Theses and Dissertations. 1325.

[2] Al-Rumaithi, Ayad, Hae-Bum Yun, and Sami F. Masri. "A Comparative Study of Mode Decomposition to Relate Next-ERA, PCA, and ICA Modes." Model Validation and Uncertainty Quantification, Volume 3. Springer, Cham, 2015. 113-133.

  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
模态分解(Modal Decomposition)是一种用于处理多模态数据的方法,它可以将复杂的数据分解成几个基本的模态分量。在Matlab中,可以使用下面的代码实现模态分解。 首先,我们需要导入相应的数据。假设我们要分解的数据存储在一个矩阵X中,其中每一列代表一个观测样本,每一行代表一个特征。 ``` % 导入数据 X = your_data; % 输入你的数据矩阵 ``` 接下来,我们可以使用奇异值分解(Singular Value Decomposition,简称SVD)来进行模态分解。SVD是一种矩阵分解的方法,可以将一个矩阵分解成三个矩阵的乘积U*S*V',其中U和V是正交矩阵,S是对角矩阵。 ``` % 变模态分解 [U,S,V] = svd(X); ``` 在SVD分解后,S矩阵的对角线元素代表了每个模态分量的能量贡献程度,可以用来判断分解的有效性。我们可以选择保留能量贡献较高的模态分量,通过设置一个能量阈值threshold,将能量贡献较低的模态分量置零。 ``` % 选择保留的模态分量 threshold = 0.95; % 设置能量阈值 total_energy = sum(diag(S).^2); % 总能量 current_energy = 0; % 当前能量 k = 0; % 保留的模态分量数 while current_energy/total_energy < threshold k = k + 1; current_energy = current_energy + S(k,k)^2; end ``` 最后,我们可以根据保留的模态分量数k,选取对应的U和V矩阵子集,并利用这些子集重构原始数据矩阵。 ``` % 重构数据 X_reconstructed = U(:,1:k) * S(1:k,1:k) * V(:,1:k)'; ``` 通过这样的步骤,我们可以将原始数据分解成k个模态分量,并根据需要进行重构。变模态分解在信号处理、图像处理等多领域有着广泛的应用。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

荔枝科研社

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值