👨🎓个人主页:研学社的博客
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。
⛳️座右铭:行百里者,半于九十。
📋📋📋本文目录如下:🎁🎁🎁
目录
💥1 概述
其中[M]、[C]、[K]分别表示系统质量矩阵、阻尼矩阵 和 刚 度 矩 阵,三 者 皆 为 实 对 称 矩 阵。{ ·· x } 、 { · x} 、{ x} 分别表示加速度、速度和位移向量。{ F ( t) } 表示系统激励力。令{ F( t) } = { F} e jωt ,{ x} = { X} e jωt ,得:
本文用于识别受高斯白噪声激励影响的2DOF系统,并增加激励和响应的不确定性(也是高斯白噪声)。
指标.EMAC : 扩展模态幅度相干
指标.MPC : 模态相位共线性
指标.CMI : 一致模式指示指标
partfac : 参与因子
📚2 运行结果
部分代码:
function [Result]=SSID(output,fs,ncols,nrows,cut)
%Input:
%output: output data of size (No. output channels, No. of data)
%fs: Sampling frequency
%ncols: The number of columns in hankel matrix (more than 2/3 of No. of data)
%nows: The number of rows in hankel matrix (more than 20 * number of modes)
%cut: cutoff value=2*no of modes
%Outputs :
%Result : A structure consist of the below components
%Parameters.NaFreq : Natural frequencies vector
%Parameters.DampRatio : Damping ratios vector
%Parameters.ModeShape : Mode shape matrix
%Indicators.EMAC : Extended Modal Amplitude Coherence
%Indicators.MPC : Modal Phase Collinearity
%Indicators.CMI : Consistent Mode Indicator
%Indicators.partfac : Participation factor
%Matrices A,C: Discrete A and C matrices
%--------------------------------------------------------------------------
[outputs,npts]=size(output); % Computes the size of matrix output
if outputs > npts % Check if output matrix size is proper or should be transposed
output=output';
[outputs,npts]=size(output);
end
%--------------------------------------------------------------------------
% Find block sizes
brows=fix(nrows/outputs); % brows = how many output blocks.
nrows=outputs*brows; % Redefine the row numbers.
bcols=fix(ncols/1); % bcols = how many time steps.
ncols=1*bcols; % Redefine the column numbers.
m=1; % inputs number.
q=outputs; % outputs number.
%--------------------------------------------------------------------------
% Form the Hankel matriices Yp,Yf.
Yp=zeros(nrows/2,ncols); %Past Output
Yf=zeros(nrows/2,ncols); %Future Output
for ii=1:brows/2
for jj=1:bcols
Yp([(ii-1)*q+1:ii*q] ,[(jj-1)*m+1:jj*m] )=output(:,((jj-1)+ii-1)*m+1:((jj)+ii-1)*m);
end
end
for ii=brows/2+1:brows
i=ii-brows/2;
for jj=1:bcols
Yf([(i-1)*q+1:i*q] ,[(jj-1)*m+1:jj*m] )=output(:,((jj-1)+ii-1)*m+1:((jj)+ii-1)*m);
end
end
%--------------------------------------------------------------------------
% Projection
O=Yf*Yp'*pinv(Yp*Yp')*Yp;
%--------------------------------------------------------------------------
% Decompose the data matrix
[R1,Sigma1,S1]=svd(O,0);
sv=diag(Sigma1);
%--------------------------------------------------------------------------
% Truncate the matrices using the cutoff
D=diag(sqrt(sv(1:cut))); % build square root of the singular values.
Dinv=inv(D); % (sigma)^(-1/2)
Rn=R1(:,1:cut); % use only the principal eigenvectors
Sn=S1(:,1:cut); % use only the principal eigenvectors
Obs=Rn*D; % Observability matrix
%--------------------------------------------------------------------------
% Calculate the realization of A and find eigenvalues and eigenvectors
A=pinv(Obs(1:nrows/2-q,:))*Obs(q+1:nrows/2,:); % build A
C=Obs(1:q,:); % build C
clear Yp Yf;
%--------------------------------------------------------------------------
% Extract the modal frequencis , damping ratios and natural frequencis
[Vectors,Values]=eig(A); % Eigenvalues and Eigenvectors
Lambda=diag(Values); % roots in the Z-plane
s=log(Lambda).*fs; % Laplace roots
zeta=-real(s)./abs(s)*100; % damping ratios
fd=(imag(s)./2./pi);
% damped natural freqs:
shapes=C*Vectors; % Mode shapes.
%--------------------------------------------------------------------------
% Calculate Modal Participation factors
InvV=inv(Vectors);
partfac=std((InvV*D*Sn')')'; %InvV involved in transforming states to the modal amplitudes by uncoupling A matrix
%--------------------------------------------------------------------------
% Calculate EMAC values
partoutput22=(C*A^(brows/2-1)*Vectors)'; %Predicted Last Block in Obs*Vectors Matrix
partoutput2=(Rn(nrows/2-q+1:nrows/2,:)*D*Vectors)'; %Actual Last Block in Obs*Vectors Matrix (Vectors involved in transforming states to the modal amplitudes by uncoupling A matrixx)
for i=1:1:cut
sum1=0;
sum2=0;
for k=1:1:q
Rik=abs(partoutput22(i,k))/abs(partoutput2(i,k));
if Rik>1
Rik=1/Rik;
end
Wik=1-(4/pi)*abs(angle(partoutput22(i,k)/partoutput2(i,k)));
if Wik<0
Wik=0;
end
EMACout=Rik*Wik;
EMACijk=EMACout;
sum1=sum1+EMACijk*(abs(partoutput2(i,k)))^2*100;
sum2=sum2+(abs(partoutput2(i,k)))^2;
end
EMACC(i)=sum1/sum2;
end
%--------------------------------------------------------------------------
% Calculate MPC values
🎉3 参考文献
部分理论来源于网络,如有侵权请联系删除。
[1]朱伟明,杨艳,刘泽远,刘程子.基于模态参与因子的无轴承开关磁阻电机振动分析与抑制[J].微电机,2022,55(10):8-15.DOI:10.15934/j.cnki.micromotors.2022.10.001.
[2]熊笑雷,赵鸣.应变模态指标在平板网架结构损伤识别分析中的应用[J].结构工程师,2012,28(02):58-65.DOI:10.15935/j.cnki.jggcs.2012.02.012.
[3] Van Overschee, Peter, and B. L. De Moor. Subspace identification for linear systems: Theory—Implementation—Applications. Springer Science & Business Media, 2012.
[4] R. Pappa, K. Elliott, and A. Schenk, “A consistent-mode indicator for the eigensystem realization algorithm,” Journal of Guidance Control and Dynamics (1993), 1993.
[5] Al Rumaithi, Ayad, "Characterization of Dynamic Structures Using Parametric and Non-parametric System Identification Methods" (2014). Electronic Theses and Dissertations. 1325.