具有模态指标的随机子空间识别【包括一致模态指标和模态参与因子】（Matlab代码实现）

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

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

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

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

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

BSRMWR 可等效为一个多自由度系统 ， 对于多自由度振动系统， 其动力学方程为

 

其中［M］、［C］、［K］分别表示系统质量矩阵、阻尼矩阵 和 刚 度 矩 阵，三 者 皆 为 实 对 称 矩 阵。{ ·· x } 、 { · x} 、{ x} 分别表示加速度、速度和位移向量。{ F ( t) } 表示系统激励力。令{ F( t) } = { F} e jωt ，{ x} = { X} e jωt ，得:

由于矩阵 ［M］ 、 ［C］ 和 ［K］ 中存在非零的非对角元素， 因此式 ( 2 ) 是一组具有耦合项的方程组 。 运用模态矩阵作为变换矩阵， 可进行坐标变换 。 其中模态矩阵和模态坐标如式( 3 ) 和 ( 4 ) 所示 。

 

本文用于识别受高斯白噪声激励影响的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.

🌈4 Matlab代码实现