System ID function AR or ARMAX
3 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
More Answers (0)
See Also
Categories
Find more on Correlation Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!