# System ID function AR or ARMAX

3 views (last 30 days)
Kin Cheong Sou on 15 Sep 2021
Commented: Kin Cheong Sou on 18 Sep 2021
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

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));
Kin Cheong Sou on 18 Sep 2021
Thank you again Ivo, for the explanation and the example. Now I understand this function much better.