ssi随机子空间算法

随机子空间算法，用于模态参数识别，可以学习

%利用随机子空间法进行结构整体模态参数识别

tic
clear
clc
% YDataz=textread（‘G:\ssi\sanxiaba5.txt‘）;%输入列向量
% YDataz=textread（‘G:\ssi\sxdq.txt‘）;%输入列向量
% YDataz=textread（‘G:\ssi\lxwlb.txt‘）;%输入列向量
YDataz=textread（‘G:\ssi\lxwlb.txt‘）;%输入列向量
% YDataz=textread（‘G:\ssi\jb4.txt‘）;%输入列向量
% YDataz=textread（‘G:\ssi\lxwdwy95lb.txt‘）;%输入列向量
YData1=YDataz‘;%YData 必须是m×n格式形式
 YData=YData1;
% YData=YData1（1:7:）;

Fs=100;
dt=1/Fs;
sampling_step=dt;
[my ny]=size（YData）;
 
%%%改变hankel总行数，即mi1的值，模态阶数保持不变，利用稳定图得到真实模态参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate hankel matrix  
% for mi1=8:2:24
mi1=42; %hankel总行数
mi2=mi1/2; %past与future分界线
mj=4096-mi1;%hankel总列数

for jj1=1:mi1   %构造块Hankel矩阵，每块是个列向量，包含了m个测点
    for jj2=1:mj
        for jj3=1:my
       Y_1_mi1（jj3+（jj1-1）*myjj2）=YData（jj3jj2+jj1-1）; 
        end
   end 
end
Yp=Y_1_mi1（1:my*mi2:）./sqrt（mj）;%分离出Yp
Yf=Y_1_mi1（（my*mi2+1）:my*mi1:）./sqrt（mj）;%分离出Yf

%对Y_1_mi1做qr分解
[Q R]=qr（Y_1_mi1‘）;
RL=R‘;%转为下三角
QT=Q‘;%Y_1_mi1值等于RL*QT
R21=RL（（my*mi2+1）:my*mi1（1:my*mi2））;%提取数据
Q1T=QT（（1:my*mi2）:）;%提取数据

% 求投影矩阵（非加权主分量算法）
OI=R21*Q1T;
% OI=Yf*Yp‘*pinv（Yp*Yp‘）*Yp;

%将OI做奇异值分解
[USV]=svd（OI）;
SS1=diag（S）;
Sum_SS1=sum（SS1）;
Length_SS1=length（SS1）;

%根据奇异熵增量来确定模态阶次 
deteE=zeros（1Length_SS1）;
E=0;
for i=1:Length_SS1
    deteE（1i）=-（SS1（i）/Sum_SS1）*log（SS1（i）/Sum_SS1）;
    E1（1i）=E+deteE（i）;
    E=E1（i）;  
end
figure;
n=10;             %模态阶次调试值
Length_Rank_S=2*n;

plot（deteE（1:Length_Rank_S））;
xlabel（‘阶次‘）;
ylabel（‘奇异熵增量‘）;
grid on
axis tight

S1=S（1:Length_Rank_S1:Length_Rank_S）;%构造奇异值非0的方阵
U1=U（1:my*mi21:Length_Rank_S）;
V1T=V（1:Length_Rank_S:）;

%可观矩阵TI
T_i=U1*（S1^0.5）;

%卡尔曼预测值，第i时刻
X_i=pinv（T_i）*OI;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%同以上原理，计算第i+1时刻卡尔曼预测值

R21_new=RL（（my*mi2+my+1）:my*mi1（1:my*mi2））;%提取数据
Q1T_new=QT（（1:my*mi2）:）;%提取数据

%求投影矩阵（非加权主分量算法）
OI_new=R21_new*Q1T_new;

%可观矩阵TI
T_i_new=T_i（1:（mi2-1）*my:）;

%卡尔曼预测值，第i+1时刻
X_i_new=pinv（T_i_new）*OI_new;
[m n]=size（X_i_new）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求系统矩阵A和C，采用最小二乘法求解
AC=[X_i_new;Y_1_mi1（（my*mi2+1）:（my*mi2+my）1:mj）]*pinv（X_i）;

WV=[X_i_new;Y_1_mi1（（my*mi2+1）:（my*mi2+my）1:mj）]-AC*X_i;

[m1 n1]=size（AC）;
A=AC（1:（m1-my）1:n1）;
C=AC（（m1-my+1）:m1:）;

%对A做特征值分解
[Eigen_Vector_A1Eigen_Value_A1]=eig（A‘nobalance‘）;
% [Eigen_Vector_A1Eigen_Value_A1]=eig（A）;
A1_diag=diag（Eigen_Value_A1）;
% 
% for iii=1:length（A1_diag）
%     lamd（iii）=log（A1_diag（iii））/dt;
%     a=real（lamd（iii））;
%     b=imag（lamd（iii））;
%     omiga（iii）=sqrt（a^2+b^2）;
%     omiga2（iii）=omiga（iii）/（2*pi）;
%     delt（iii）=-a/sqrt（a^2+b^2）;
% end

% F1=omiga2;
% D1=delt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求模态参数
% 计算模态频率向量
F1=abs（log（A1_diag‘））./（2*pi*dt）;

% 计算阻尼比向量
D1=sqrt（1./（（（imag（log（ A1_diag‘））./real（log（ A1_diag

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3758  2010-06-09 12:54  ssi.m

----------- ---------  ---------- -----  ----

                 3758                    1