ssi—1.zip_SSI法_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)*my,jj2)=YData(jj3,jj2+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做奇异值分解 [U,S,V]=svd(OI); SS1=diag(S); Sum_SS1=sum(SS1); Length_SS1=length(SS1); %根据奇异熵增量来确定模态阶次 deteE=zeros(1,Length_SS1); E=0; for i=1:Length_SS1 deteE(1,i)=-(SS1(i)/Sum_SS1)*log(SS1(i)/Sum_SS1); E1(1,i)=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_S,1:Length_Rank_S);%构造奇异值非0的方阵 U1=U(1:my*mi2,1: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_A1,Eigen_Value_A1]=eig(A,'nobalance'); % [Eigen_Vector_A1,Eigen_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'))).^2)+1)); %计算振型 ZX1=C*Eigen_Vector_A1; ZX2=imag(ZX1);%取虚部，但有重根 % ZX=ZX2(1:2:end); [F2,I]=sort(F1); %剔除方程解中的非模态项(非共轭根)和共轭项(重复项) mm=0; for k=1:(length(A1_diag)-1) if F2(k)~=F2(k+1) continue; end mm=mm+1; l=I(k)'; F(mm)=F1(l); %自振频率 D(mm)=D1(l);%阻尼比 end % end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% toc