【参数辨识】随机子空间辨识 (SSI)(Matlab实现)

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

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

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

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现


💥1 概述

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

基本原理:

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

  1. 数据预处理

    • 收集系统的输入和输出数据。
    • 对数据进行预处理,如去除均值、归一化等,以提高后续计算的准确性。
  2. 构造观测矩阵

    • 从输入-输出数据中构造两个观测矩阵 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)

  3. 奇异值分解 (SVD)

    • 对观测矩阵 Y1Y1 进行奇异值分解,得到:

      Y1=UΣVTY1=UΣVT

    • 选择前 nn 个最大的奇异值对应的左奇异向量组成矩阵 UnUn,这些向量构成了系统的状态空间基。
  4. 系统矩阵估计

    • 利用 UnUn 和 Y2Y2 计算系统矩阵 AA 和 CC:

      A^=Un+Y2VnΣn−1A^=Un+Y2VnΣn−1

      C^=Y1VnΣn−1C^=Y1VnΣn−1

    • 其中 Un+Un+ 表示 UnUn 的伪逆。
  5. 模态参数提取

    • 通过对系统矩阵 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.

🌈4 Matlab代码实现

图片

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值