System ID function AR or ARMAX

3 views (last 30 days)
Hello,
I am trying to fit a sum of two exponentials with time-domain data using system identification toolbox functions ar or armax. It appears that the fitting quality is not good and it's not my expectation. I attach my simple code below. I wonder if I made any simple mistakes. Any idea would be welcome. Thanks!
% generate data: sum of 2 exp
N = 1000;
x = (0:(N-1))';
y0 = exp(-0.01*x)+0.8*exp(-0.1*x); % positive sum of two decaying exp
y = y0+0.01*max(y0)*randn(N,1); % add a bit of noise
% apply SYSID tools: ar or armax
n = 2;
yid = iddata(y,[],1);
%sys = ar(yid,n);
sys = armax(yid,[n 0]);
[yh, FIT, x0] = compare(yid,sys);
% show results
figure(1);
clf;
plot(x,y,x,yh.y,'LineWidth',4);
legend('data','SYSID');
disp('system poles');
disp(roots(sys.a));
Kin

Accepted Answer

Ivo Houtzager
Ivo Houtzager on 16 Sep 2021
Edited: Ivo Houtzager on 17 Sep 2021
Define a step input, and fit OE model as this model structure fits your data the best. See updated code below.
% generate data: sum of 2 exp
d = 10; % input delay
p = 100; % start step
N = 1000;
x = [zeros(p+d,1); (1:(N-p-d))'];
y0 = exp(-0.01*x)+0.8*exp(-0.1*x); % positive sum of two decaying exp
y = y0+0.01*max(y0)*randn(N,1); % add a bit of noise
y = y - y0(1); % remove start value, as fit model around steady state zero
u = [zeros(p,1); ones(N-p,1)];
yid = iddata(y,u,1);
% apply SYSID tools: oe
n = 2;
[sys1,ic1] = oe(yid,[n n d]);
[yh1, FIT1] = compare(yid,sys1);
% apply SYSID tools: oe with initial estimate
init_sys = idpoly(zpk(0.8,[0.99 0.9],-y0(1),1)); % initial sys
init_sys.nk = d; % number of samples input delay
opt = oeOptions('InitialCondition','zero'); % initial state
[sys2,ic2] = oe(yid,init_sys,opt);
[yh2, FIT2] = compare(yid,sys2);
% show results
figure(1);
clf;
t = 0:N-1;
plot(t,u,t,y+y0(1),t,yh1.y+y0(1),t,yh2.y+y0(1),'LineWidth',4);
legend('input','data','SYSID1','SYSID2');
disp('system poles');
disp(roots(sys2.f));
  3 Comments
Kin Cheong Sou
Kin Cheong Sou on 18 Sep 2021
Thank you again Ivo, for the explanation and the example. Now I understand this function much better.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!