ERA辨识的问题

quaintchewster

我不是学控制的，所以对这些玩意是新手。最近做有关ERA辨识，今天弄了一个小玩意测试一下。

ERA 代码是网上下的，是SISO的。我自己用一个系统生成了step响应数据，调用era进行辨识。得到的系统画出的step响应和impulse响应都和原系统很吻合。但是矩阵和传递函数看起来很不一样啊。矩阵本不是唯一，但是传递函数是不是唯一的？有什么别的方法可以进一步验证呢？

%-- 8/5/07 1:43 PM --%
k1=10;
k2=50;
b1=1;b2=2;
m1=1;m2=1;
A=[0 0 1 0;0 0 0 1;-k1/m1 k1/m1 -b1/m1 b1/m1;k1/m2 -(k1+k2)/m2 b1/m2 -(b1+b2)/m2]
B=[0 0 1/m1 0]';
C=[1 0 0 0];
D=[0];
[ys,x,t]=step(A,B,C,D);
yp=ys;
% 由step响应构造impulse响应
for i=1:size(ys)
if(i>1), yp(i) = ys(i)-ys(i-1);end
end

% yp 信号
% 4 阶数
% 200 信号长度
% t(2)-t(1) 信号时间步长
% 2 要求输出连续系统矩阵
[cA,cB,cC,cD]=era(yp,4,200,t(2)-t(1),2);
cA
step(cA,cB,cC,cD)

% 传递函数
[onum,oden]=ss2tf(A,B,C,D)
[cnum,cden]=ss2tf(cA,cB,cC,cD)

% step响应
[cy,cx,ct]=step(cA,cB,cC,cD);
[cy,cx,ct]=step(cA,cB,cC,cD);
[oy,ox,ot]=step(A,B,C,D);
pstep=plot(ot,oy,'*',ct,cy,'-')

% impulse响应
[cy,cx,ct]=impulse(cA,cB,cC,cD);
[oy,ox,ot]=impulse(A,B,C,D);
pimp=plot(ot,oy,'*',ct,cy,'-')

function [A,B,C,D]=era(h,n,N,Ts,def);

% Eigensystem Realization Algorithm (ERA)
%
% Author: Samuel da Silva - UNICAMP
% e-mail: samsilva@fem.unicamp.br
% Date: 2006/10/20
%
% [A,B,C,D]=era(h,n,N,Ts,def);
%
% Inputs:
% h: discrete-time impulse response
% n: order of the system
% N: number of samples to assembly the Hankel matrix
% Ts: sample time
% def: if = 1: the output will be the discrete-time state-space model
% if = 2: the output will be the continuous-time state-space model
%
% Otputs:
% [A,B,C,D]: state-space model
%
% Note: For now, it works to SISO systems and it is necessary the control toolbox
%
% References: Juang, J. N. and Phan, M. Q. "Identification and Control of
% Mechanical Systems", Cambridge University Press, 2001

% Hankel matrix
H0 = hankel(h(2:N+1)); % k = 0
H1 = hankel(h(3:N+2)); % k = 1;

% Factorization of the Hankel matrix by use of SVD
[R,Sigma,S] = svd(H0);
% R and S are orthonormal and Sigma is a rectangular matrix

Sigman = Sigma(1:n,1:n);

Wo = R(:,1:n)*Sigman^0.5; % observability matrix
Co = Sigman^.5*S(:,1:n)'; % controllability matrix

% The identified system matrix are:
A = Sigman^-.5*R(:,1:n)'*H1*S(:,1:n)*Sigman^-.5; % dynamic matrix
B = Co(:,1); % input matrix
C = Wo(1,:); % output matrix
D = h(1); % direct-transmission matrix

sysdisc = ss(A,B,C,D,Ts); % discrete-time system

if def == 2
syscont = d2c(sysdisc,'zoh'); % Conversion of discrete LTI models to continuous time
[A,B,C,D]=ssdata(syscont); % continuous system
end

%--------------------------------------------------------------------------