System ID function AR or ARMAX

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

Answer 1

[Ivo Houtzager]

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]

Thanks Ivo. The fitting result looks much better now. Though I encounter, once in a while, some unexpected results like

system poles

0.9873

-0.5199

where the expected ones should be exp(-0.01) \approx 0.99 and exp(-0.1) \approx 0.9048. I will look deeper into this.

Another issue: I do get some warning (below) while running your provided code. I try to pad some zeros before the ones in the input u and pad some 1.8's (or 0's) before y, but the results are not expected. I wonder if you have any idea how to change the code to avoid the warning?

Warning: For transient data (step or impulse experiment), make sure that the change in input signal

does not happen too early relative to the order of the desired model. You can achieve this by

prepending sufficient number of zeros (equilibrium values) to the input and output signals. For

example, a step input must be represented as [zeros(nx,1); ones(N,1)] rather than ones(N,1), such

that nx > model order.

[Ivo Houtzager]

Compared to the ARX model, the OE model optimization problem is non-linear. Thus the optimization is non-convex. Depending on the initial values, the non-linear optimization can get sometimes stuck in local optimum instead of going to global optimum. If you have an initial estimate or guess of parameters, you can give that as an option to the OE function to improve the convergence. The ARX model optimization problem is linear and does not have this issue as this optimization is convex, but will require high-order model to get a decent fit on your data.

For your second issue, you also have to pad the same number of samples of constant values to the output y as well.

I have updated the code in answer with example to remove the warning and to apply an initial estimate.

[Kin Cheong Sou]

Thank you again Ivo, for the explanation and the example. Now I understand this function much better.