How to get good sys ids using n4sid

5 views (last 30 days)
Alan
Alan on 15 Nov 2015
I have a test problem to try out the n4sid function in the System Identification Toolbox.
I have simulated a simple SISO second order system. The only interesting feature is that it is proper but not strictly proper. I can get a good fit in reasonable time using tfest, but I'd like to use n4sid to learn how to use it and get good results. (ssest takes too long on this problem.)
The code below produces this graph, with test taking 28 seconds and n4sid taking 32 seconds.
ftest is almost spot on, but as you can see, n4sid is not close to the correct solution.
What should I do to tweak n4sid to get a better solution? A secondary question is what do I do to make n4sid return a continuous time model? The help says "For estimating a continuous-time model, specify 'Ts'/0 as name-value pair" which is what I thought I did, but the model returned is discrete time.
% generate a simple SISO SOS
% this is not strictly proper, but is proper
% example, series RCL circuit, L=1 H, R=0.1 ohm, C=1 F
% transfer function from input voltage to voltaga across inductor
m = idtf([1 0 0],[1 0.1 1]);
% generate a random input signal
N=1e6;
u=randn(N,1);
[B,A] = butter(4,0.75);
u=filter(B,A,u);
% assume a 100 Hz sample rate
ts=.01;
fs=1/ts;
w0=linspace(0,pi*fs,5001);
% simulate the output
z=iddata([],u,ts);
y=sim(m,z);
% add some noise to the output
e=0.02*randn(N,1);
e=filter(B,A,e);
y=y.OutputData+e;
% create iddata structure
z=iddata(y,u,ts);
% tf model using test
% fix the structure to ensure double zero at origin
init_sys = idtf([NaN 0 0],[1 NaN NaN]);
init_sys.Structure.num.Free([2,3]) = false;
disp('tfest2')
tic
m4=tfest(z,init_sys);
T4=freqresp(m4,w0);
T4=squeeze(T4);
toc
% get ss model using n4sid
disp('n4sid')
tic
m6=n4sid(z,4,'Feedthrough',true,'Form','Companion','Ts',0);
toc
T6=freqresp(m6,w0);
T6=squeeze(T6);
% plot all
loglog(w0,abs(T4),w0,abs(T6));
legend('tfest','n4sid')
grid on
figure(gcf())

Answers (0)

Community Treasure Hunt

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

Start Hunting!