quaintchewster 发表于 2007-8-5 15:54

ERA辨识的问题

我不是学控制的,所以对这些玩意是新手。最近做有关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=
B=';
C=;
D=;
=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 要求输出连续系统矩阵
=era(yp,4,200,t(2)-t(1),2);
cA
step(cA,cB,cC,cD)

% 传递函数
=ss2tf(A,B,C,D)
=ss2tf(cA,cB,cC,cD)

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

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

function =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
%
% =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:
% : 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
= 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
=ssdata(syscont); % continuous system
end

%--------------------------------------------------------------------------
页: [1]
查看完整版本: ERA辨识的问题