应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 专业学科区 » 信息科学 » 控制工程 » MIMO系统子空间系统辨识算法求助
101/1返回列表
查看: 2464  |  回复: 9
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

书寻玉

铁杆木虫 (正式写手)

[求助] MIMO系统子空间系统辨识算法求助已有1人参与

我想问下，大家有没有做过MIMO系统的子空间系统辨识啊，最好是双输入双输出系统的辨识，有没有相关的资料共享下，最好是MATLAB的例程。我按书上的方法编制了MATLAB程序，算出来不对，其中B、D怎么用最小二乘法求呢，最好有详细的资料提供下啊，不胜感激啊。
回复此楼

» 猜你喜欢
1楼2014-01-03 10:36:41
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

southisland

新虫 (初入文坛)

【答案】应助回帖

我写的一段程序，可以参考下
/////////////////////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
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

southisland

新虫 (初入文坛)

【答案】应助回帖

找到了当时写的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+12*m+p)*I,m*I+12*m+p)*I))...
        *R(m*I+12*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(1I-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
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

书寻玉

铁杆木虫 (正式写手)
引用回帖:
2楼: Originally posted by southisland at 2014-02-23 16:50:39
我写的一段程序，可以参考下
/////////////////////step 5 利用最小二乘计算B、D///////////////////
                /*参考Subspace Identification for Linear Systems:
                Theory Implementation Applications,pp.56*/
                M ...

非常感谢你，能分享下其他的step吗？非常感谢啊
赞一下 回复此楼
4楼2014-02-23 18:48:04
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

southisland

新虫 (初入文坛)

【答案】应助回帖

引用回帖:
4楼: Originally posted by 书寻玉 at 2014-02-23 18:48:04
非常感谢你，能分享下其他的step吗？非常感谢啊...

没有其他步骤了，基本的SIM算法就这几个步骤的。

[ 发自小木虫客户端 ]
赞一下 回复此楼
5楼2014-02-24 16:45:14
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

书寻玉

铁杆木虫 (正式写手)
引用回帖:
5楼: Originally posted by southisland at 2014-02-24 16:45:14
没有其他步骤了，基本的SIM算法就这几个步骤的。
...

嗯，谢谢了，是SIMO的吗？不是MIMO的吗
赞一下 回复此楼
6楼2014-02-24 21:33:50
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

southisland

新虫 (初入文坛)

【答案】应助回帖

引用回帖:
6楼: Originally posted by 书寻玉 at 2014-02-24 21:33:50
嗯，谢谢了，是SIMO的吗？不是MIMO的吗...

子空间辨识算法最大优点在于不用优化求解，而且处理多变量系统的辨识时算法无需修改，只是将辨识数据组成向量即可，是一种真正意义上的多变量系统辨识算法。

[ 发自小木虫客户端 ]
赞一下 回复此楼
7楼2014-02-24 21:50:24
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

秦时小沫

新虫 (小有名气)
引用回帖:
2楼: Originally posted by southisland at 2014-02-23 16:50:39
我写的一段程序，可以参考下
/////////////////////step 5 利用最小二乘计算B、D///////////////////
                /*参考Subspace Identification for Linear Systems:
                Theory Implementation Applications,pp.56*/
                M ...

您好，我想问一下，如何将闭环子空间辨识算法应用到串级控制系统啊，此时我们如果要辨识内环的模型，如何辨识呢，望您指导，感激不尽！
赞一下 回复此楼
8楼2014-03-17 14:54:55
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

southisland

新虫 (初入文坛)

【答案】应助回帖

引用回帖:
8楼: Originally posted by 秦时小沫 at 2014-03-17 14:54:55
您好，我想问一下，如何将闭环子空间辨识算法应用到串级控制系统啊，此时我们如果要辨识内环的模型，如何辨识呢，望您指导，感激不尽！...

从算法上看，闭环SIM辨识需要三组数据，即参考输入信号r，系统输出信号y和控制器输出信号u，基于这些测量得到的值可以使用闭环算法得到内环模型。具体的实验设计可以参考相关论文和书籍。

[ 发自小木虫客户端 ]
赞一下 回复此楼
9楼2014-03-21 09:21:33
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lovskyboy416

新虫 (初入文坛)
您好，您的mimo子空间辨识算法能发给一下做个参考吗

发自小木虫Android客户端
赞一下 回复此楼
10楼2017-05-30 20:32:36
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 书寻玉 的主题更新
101/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭