| 查看: 2594 | 回复: 9 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
书寻玉铁杆木虫 (正式写手)
|
[求助]
MIMO系统子空间系统辨识算法求助 已有1人参与
|
|
| 我想问下，大家有没有做过MIMO系统的子空间系统辨识啊，最好是双输入双输出系统的辨识，有没有相关的资料共享下，最好是MATLAB的例程。我按书上的方法编制了MATLAB程序，算出来不对，其中B、D怎么用最小二乘法求呢，最好有详细的资料提供下啊，不胜感激啊。 |
回复此楼» 猜你喜欢
生物学308求调剂（一志愿华东师大）
已经有6人回复
-
调剂求收留
已经有7人回复
-
272分材料子求调剂
已经有36人回复
-
275求调剂
已经有8人回复
-
一志愿211，化学学硕，310分，本科重点双非，求调剂
已经有20人回复
-
070300化学学硕311分求调剂
已经有19人回复
-
材料与化工调剂
已经有13人回复
-
材料与化工调剂
已经有33人回复
-
复试调剂
已经有7人回复
-
一志愿哈工大 085600 277 12材科基求调剂
已经有17人回复
书寻玉
铁杆木虫 (正式写手)
- 应助: 2 (幼儿园)
- 金币: 5543.2
- 散金: 3
- 红花: 22
- 帖子: 478
- 在线: 241.1小时
- 虫号: 2209924
- 注册: 2012-12-27
- 性别: GG
- 专业: 制造系统与自动化
4楼2014-02-23 18:48:04
【答案】应助回帖
|
我写的一段程序，可以参考下 /////////////////////step 5 利用最小二乘计算B、D/////////////////// /*参考[Peter Van Overschee,96]Subspace Identification for Linear Systems: Theory Implementation Applications,pp.56*/ Matrix<Type> Mi,Li,MM(I*(p*I-n),m),LL(I*(p*I-n),p*I),IGamd(p*I,p+n),BD; Mi=trT(U2)*R31*pinv(R11); Li=trT(U2); for (int i=1;i<=I;i++) { for(int j=1;j<=p*I-n;j++) { for(int k=1;k<=m;k++) MM((i-1)*(p*I-n)+j,k)=Mi(j,(i-1)*m+k); } } for(int i=1;i<=I*(p*I-n);i++) for(int j=1;j<=p*I;j++) LL(i,j)=double(0); for (int i=1;i<=I;i++) { for(int j=i;j<=I;j++) { for(int k=1;k<=p*I-n;k++) { for(int t=1;t<=p;t++) LL((i-1)*(p*I-n)+k,(j-i)*p+t)=Li(k,(j-1)*p+t); } } } //设置矩阵IGamd=[Ip,0;0,Gamd] for (int i=1;i<=p;i++) { for(int j=1;j<=p+n;j++) { if (j==i) IGamd(i,j)=double(1); else IGamd(i,j)=double(0); } } for (int i=1;i<=p*(I-1);i++) { for(int j=1;j<=p;j++) IGamd(p+i,j)=double(0); for(int j=1;j<=n;j++) IGamd(p+i,p+j)=Gamd(i,j); } //QRD<Type> qrBD;qrBD.dec(LL*IGamd);BD=qrBD.solve(MM); BD=pinv(LL*IGamd)*MM; //最小二乘计算BD=[D;B] for(int i=0;i<n+p;i++) { if(i<p) D.setRow(BD.getRow(i),i); //提取D矩阵 else B.setRow(BD.getRow(i),i-p); //提取B矩阵 } //计算算法耗时 |
2楼2014-02-23 16:50:39
【答案】应助回帖
|
找到了当时写的matlab函数 function [A,B,C,D]=mysubid(Y,U,I,obqr) %基本的线性开环子空间辨识算法 %Y：辨识输出信号， U：辨识输入信号， I：Hankel矩阵行数 %obqr:如果使用QR分解的斜向投影，则置于非0的任意数 %obqr==0,使用斜向投影 %计算时可选直接投影或者基于QR分解的斜向投影 %extract matrices A & C from matric Gam,then Solve least square: %M=L*X*[B;D] % 2012-4-11 rewrite 2012-5-11 %2012-7-11 增加可测噪声模型的辨识 %%%START%%% if (nargin < 4);obqr = 0;end [m,N]=size(U);p=size(Y,1); J=N-2*I+1; % obqr = 1; %% STEP 1: 构造输入输出数据hankel矩阵块% HklU=blkhank(U,I,J ); Up=[];Yp=[];Uf=[];Yf=[]; for i=1:I Up=[Up;U(:,i:J+i-1)]; Uf=[Uf;U(:,I+i:J+I+i-1)]; Yp=[Yp;Y(:,i:J+i-1)]; Yf=[Yf;Y(:,I+i:J+I+i-1)]; end %% STEP 2：投影并加权（直接投影和QR分解斜向投影） %直接正交投影 Wp=[Yp;Up]; if obqr==0 PI_Uf=eye(J)-Uf'*pinv(Uf*Uf')*Uf; O=Yf*PI_Uf*Wp'*pinv(Wp*PI_Uf*Wp')*Wp; else %采用QR分解计算斜向投影 YuW=[Uf;Wp;Yf]; [yu_q,yu_r]=qr(YuW'); R= yu_r';% R = R(1:2*I*(m+p),1:2*I*(m+p)); % Truncate Q=yu_q'; O=R(I*p+2*I*m+1:2*I*(m+p),p*I+1:I*(2*p+m))*inv(R(m*I+1  2*m+p)*I,m*I+12*m+p)*I))...*R(m*I+1 2*m+p)*I,12*m+p)*I)*Q(12*m+p)*I, ;end %% STEP 3：SVD分解，选择系统状态阶数 [Uu,S,V] = svd(O); ss = diag(S)'; lnss=log(ss); h = bar(1:I*p,lnss,'b'); grid on;title('Singular Value Histogram');xlabel('The Number of Singular Value ');ylabel('Singular Value'); n = input(' System order ? ');%选择系统阶数 %% STEP 4：由Gam利用其移不变性计算C和A U1=Uu(:,1:n);V1=V(1:n, ;r=sqrt(S(1:n,1:n)); Gam=U1*r; Xf=r*V1; C_h=Gam(1:p,1:n); A_h=pinv(Gam(1:p*(I-1),1:n))*Gam(p+1:p*I,1:n); %% STEP 5：计算B和D: 对于M=L*X*[B;D]采用LS计算B,D tao=Gam';%注意：转置矩阵 PI_Gam=eye(size(tao,2))-tao'*inv(tao*tao')*tao;%Gam转置的正交投影算子，PI_Gam*Gam=0 Ms=PI_Gam*Yf*pinv(Uf); for k=1:I M(I*p*(k-1)+1:I*p*k, =Ms(:,m*(k-1)+1:m*k);%构造M的列块矩阵L(I*p*(k-1)+1:I*p*k, =[PI_Gam(:,p*(k-1)+1:p*I),zeros(p*I,p*(k-1))];%构造L的块方阵end IG = [eye(p),zeros(p,n);zeros(p*(I-1),p),Gam(1 I-1)*p,];%构造【I,0;0,gamma】矩阵% Solve least squares sol_bd = (L*IG)\M; D_h = sol_bd(1:p, ;B_h = sol_bd(p+1:p+n,; |
3楼2014-02-23 16:55:40
5楼2014-02-24 16:45:14
