rels.m

系统辨识的M文件

function thita=rels(z, na, nb, nc)
% 用增广矩阵法估计参数 ....                                      % 函数说明体
if isa(z,'iddata')                                               % 输入输出数据
   y = pvget(z,'OutputData');
   u = pvget(z,'InputData');
   z = [y{1},u{1}];
end
[nz,ns]=size(z);     nn=na+nb+nc;                                % 数据维数与阶次
thitak=ones(nn,1)*0.001;                                         % 给定初始条件
thita=zeros(nn,nz);   
p1=eye(nn,nn)*(1.0e6);    p2=zeros(nn,nn);
K=zeros(nn,1);           e=zeros(nz,1);
I=eye(nn,nn);
for i=na+1:nz
    Q=[[-z(i-1:-1:i-na, 1)]',[z(i-1:-1:i-nb, 2)]',[e(i-1:-1:i-nc, 1)]']; % 数据向量   
    K=p1*Q'/(Q*p1*Q'+1);                                         % 增益因子
    p2=(I-K*Q)*p1;                                               % 协方差阵
    thita(:,i)=thitak+K*(z(i,1)-Q*thitak);                       % 参数估计值
    p1=p2;
    thitak=thita(:,i);
    e(i)=z(i,1)-Q*thitak;                                        % 预报误差
end