SSI.zip_SSI_SSI matlab_SSI程序_ssi matlab_改进SVD

%nsys=3; % nsys=系统自由度 %Mss=eye(3)*4*10^5;% Mss=结构质量矩阵(kg) %Kss=[4,-2,0;-2,4,-2;0,-2,2]*10^8;% Kss=结构刚度矩阵(N/m) %css=[1.3506,-0.5286,0;-0.5286,1.3506,-0.5286; % 0,-0.5286,0.8220]*10^4;%css=结构阻尼矩阵(N.s/m) nsys=3; % nsys=系统自由度 Mss=eye(3);% Mss=结构质量矩阵 Kss=[2,-1,0;-1,2,-1;0,-1,1];% Kss=结构刚度矩阵 css=[0.151,-0.0591,0;-0.0591,0.151,-0.0591; 0,-0.0591,0.0919];%css=结构阻尼矩阵 dt=0.5;%采样时间 Ca=eye(3);Cv=zeros(3,3);Cd=zeros(3,3); % 输出=Ca*加速度+Cv*速度+Cd*位移，这里只输出加速度 Ass = [ zeros(nsys,nsys),eye(nsys,nsys); -inv(Mss)*Kss,-inv(Mss)*css]; % A为状态矩阵，表示系统内部状态的联系 B2=zeros(nsys,1);B2(1,1)=1; % B2为描述输入位置的矩阵，设置为在第一个自由度节点处激励 Bss=[zeros(nsys,nsys);inv(Mss)]*B2; % B为输入矩阵，表示输入对系统的作用关系 Css = [Cd-Ca*Mss\Kss,Cv-Ca*Mss\css]; % C为输出矩阵，表示输出与系统内部状态之间的关系 Dss=Ca*Mss\B2; % D为传递系数矩阵，表示输入对输出的作用关系 bareStruct=ss(Ass,Bss,Css,Dss); t=0:dt:1000; [Yt,T]=impulse(bareStruct,t); Y=Yt' %************************************************************* %HY为输出Hankel矩阵 %R用于表示Hankel矩阵的块行数目，是需人为确定的一个系数。它至少要大于想要 %计算的最大系统计算阶次。理论上仅需要大于能观性指数，但由于能观性指 %数是未知的，通常假设大于系统的阶次。但的选取不能选的太大,因为计算 %机计算时间的长短与R的平方成正比。 R=4; [nr, nc] = size(Y); % Y的列为同一时刻各个传感器的输出，Y的行为同一传感器各个时刻的输出 HY = zeros(nr*R, nc-R+1); for r=1:R, for c=1:nc-R+1 HY(1+nr*(r-1):nr*r,c) = Y(:,r+c-1); end end HY; %************************************************************* %HYp为输出Hankel矩阵的上半块 [nr, nc] = size(Y); [HYr,HYc]=size(HY); HYp=zeros(HYr/2,nc-R+1); HYp=HY(1:HYr/2,:); HYp; %************************************************************* %HYf为输出Hankel矩阵的下半块 [nr, nc] = size(Y); [HYr,HYc]=size(HY); HYf=zeros(HYr/2,nc-R+1); HYf=HY(HYr/2+1:HYr,:); HYf; %************************************************************* %HYp_p为输出Hankel矩阵的上半块加上下半块的第一行块 [nr, nc] = size(Y); [HYr,HYc]=size(HY); HYp_p=zeros(HYr/2+nr,nc-R+1); HYp_p=HY(1:HYr/2+nr,:); HYp_p; %************************************************************* %HYf_m为输出Hankel矩阵的下半块减去下半块的第一行块 [nr, nc] = size(Y); [HYr,HYc]=size(HY); HYf_m=zeros(HYr/2-nr,nc-R+1); HYf_m=HY(HYr/2+nr+1:HYr,:); HYf_m; %************************************************************* %这里U采用Matlab-impulse输入的冲击信号，单点激励 [nr, nc] = size(Y); U=zeros(1,nc);U(1,1)=1; %HU为输入Hankel矩阵 % R为Hankel矩阵的行块数,R取偶数 [ur, uc] = size(U); %U的列为同一时刻各个激励点的输入，U的行为同一激励点各个时刻的输入 HU = zeros(ur*R, uc-R+1); for r=1:R, for c=1:uc-R+1 HU(1+ur*(r-1):ur*r,c) = U(:,r+c-1); end end HU; %************************************************************* %HUp为输入Hankel矩阵的上半块 [ur, uc] = size(U); [HUr,HYc]=size(HU); HUp=zeros(HUr/2,uc-R+1); HUp=HU(1:HUr/2,:); HUp; %************************************************************* %HUf为输入Hankel矩阵的下半块 [ur, uc] = size(U); [HUr,HYc]=size(HU); HUf=zeros(HUr/2,uc-R+1); HUf=HU(HUr/2+1:HUr,:); HUf; %************************************************************* %HUp_p为输入Hankel矩阵的上半块加上下半块第一行块矩阵 [ur, uc] = size(U); [HUr,HYc]=size(HU); HUp_p=zeros(HUr/2+ur,uc-R+1); HUp_p=HU(1:HUr/2+ur,:); HUp_p; %************************************************************* %HUf_m为输入Hankel矩阵的下半块减去下半块的第一行块 [ur, uc] = size(U); [HUr,HYc]=size(HU); HUf_m=zeros(HUr/2-ur,uc-R+1); HUf_m=HU(HUr/2+ur+1:HUr,:); HUf_m; %************************************************************* %HWp上块为HUp矩阵，下块为HYp矩阵 [HUpr,HUpc]=size(HUp); [HYpr,HYpc]=size(HYp); HWp=zeros(HUpr+HYpr,HUpc); HWp(1:HUpr,:)=HUp(:,:); HWp(HUpr+1:HUpr+HYpr,:)=HYp(:,:); HWp; %************************************************************* %HWp_p上块为HUp_p矩阵，下块为HYp_p矩阵 [HUp_pr,HUp_pc]=size(HUp_p); [HYp_pr,HYp_pc]=size(HYp_p); HWp_p=zeros(HUp_pr+HYp_pr,HUp_pc); HWp_p(1:HUp_pr,:)=HUp_p(:,:); HWp_p(HUp_pr+1:HUp_pr+HYp_pr,:)=HYp_p(:,:); HWp_p; %************************************************************* %************************************************************* %************************************************************* %%%求gam_i [HWpr,HWpc]=size(HWp); [HYfr,HYfc]=size(HYf); [Before_pseudo_inverse_r,Before_pseudo_inverse_c]... =size([HWp*HWp.',HWp*HUf.';HUf*HWp.',HUf*HUf.']); %O_i为计算HYf通过HUf在HWp上的斜投影，pinv为求广义逆矩阵， O_i=HYf*[HWp.',HUf.']*pinv([HWp*HWp.',HWp*HUf.';HUf*HWp.',HUf*HUf.'])... *[eye(HWpr);zeros(Before_pseudo_inverse_c-HWpr,HWpr)]*HWp; %对O_i进行加权，这里使用最简单的UPC算法进行加权,然后进行SVD W1=eye(HYfr);W2=eye(HYfc); WO_iW=W1*O_i*W2; [U,S,V]=svd(WO_iW); %从U与S中提取相应的U1和S1 index=find(S>0.00001); S1=diag(S(index)); [S1r,S1c]=size(S1); U1=U(:,1:S1c); gam_i=inv(W1)*U1*sqrt(S1); %************************************************************* %************************************************************* %************************************************************* %%%求gam_i_1 [HWp_pr,HWp_pc]=size(HWp_p); [HYf_mr,HYf_mc]=size(HYf_m); [Before_pseudo_inverse_1_r,Before_pseudo_invers_1_c]... =size([HWp_p*HWp_p.',HWp_p*HUf_m.';HUf_m*HWp_p.',HUf_m*HUf_m.']); %O_i_1为计算HYf_m通过HUf_m在HWp_p上的斜投影，pinv为求广义逆矩阵， O_i_1=HYf_m*[HWp_p.',HUf_m.']... *pinv([HWp_p*HWp_p.',HWp_p*HUf_m.';HUf_m*HWp_p.',HUf_m*HUf_m.'])... *[eye(HWp_pr);zeros(Before_pseudo_invers_1_c-HWp_pr,HWp_pr)]*HWp_p; %对O_i_1进行加权，这里使用最简单的UPC算法进行加权,然后进行SVD W1_1=eye(HYf_mr);W2_1=eye(HYf_mc); W_1O_i_1W_1=W1_1*O_i_1*W2_1; [U_1,S_1,V_1]=svd(W_1O_i_1W_1); %从U_1与S_1中提取相应的U1_1和S1_1 index=find(S_1>0.00001); S1_1=diag(S_1(index)); [S1_1r,S1_1c]=size(S1_1); U1_1=U_1(:,1:S1_1c); gam_i_1=inv(W1_1)*U1_1*sqrt(S1_1); %************************************************************* %计算Xi和Xip Xi=pinv(gam_i)*O_i; Xip=pinv(gam_i_1)*O_i_1; %************************************************************* %%%求得A，B，C，D,模型建立值为Ass,Bss,Css,Dss %[Xip;HY(HYr/2+1:HYr/2+nr,:)]=[A,B;C,D]*[Xi;HU(HUr/2+1,:)] L=[Xip;HY(HYr/2+1:HYr/2+nr,:)]; R=[Xi;HU(HUr/2+1,:)]; sol=(pinv(R*R')*R*L')'; [solr,solc]=size(sol); A=sol(1:S1r,1:S1c); B=sol(1:S1r,S1c+1:solc); C=sol(S1r+1:solr,1:S1c); D=sol(S1r+1:solr,S1c+1:solc); %真实A、B、C、D %Ass,Bss,Css,Dss %************************************************************* %真实固有频率、阻尼比 ex=Kss*inv(Mss); [ass, lamss] = eig(ex); [lamssr,lamssc]=size(lamss); lamss_r=lamss(1:lamssr,1:lamssc); realss_r=real(lamss_r); imagss_r=imag(lamss_r); wrss=sqrt(realss_r^2+imagss_r^2) kexiss=abs(realss_r)/wrss %************************************************************* %识别固有频率、阻尼比 [a, lam] = eig(A); [lamr,lamc]=size(lam); lam lam_r=diag(lam(1:2:lamr,1:2:lamc)) lam_rc=log(lam_r)/dt real_r=real(lam_rc); imag_r=imag(lam_rc); wr=sqrt(real_r.^2+imag_r.^2); wr=diag(wr) kexi=abs(diag(real_r))/wr %************************************************************* %************************************************************* rank(O_i)