%利用隨機子空間法進行結構整體模態參數識別

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