Grid search program to determine system transfer function

determineSystem.m

function besttf = determineSystem(z, npole, nzero)
%% determineSystem.m
% Grid search a system model from the given number of zeros and poles
%
%% Inputs: 
%   z      : iddata object (created by iddata(y,u,Ts))
%   npole  : maximum number of poles for grid search
%   nzero  : maximum number of zeros for grid search
% 
%% Outputs:
%   besttf : Transfer function which is most likely 
% 
%% Example: 
% (Pitch-elevator transfer function:)
% (http://ctms.engin.umich.edu/CTMS/index.php?example=AircraftPitch&section=SystemModeling)
% s = tf('s');
% G = (1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s); 
% ts=0.1; t=0:ts:20; nx=10;
% u=[zeros(nx,1)',ones(length(t)-nx,1)'];
% [y,t]=lsim(G,u,t);
% z=iddata(y,u',0.1);
% determineSystem(z, 6, 6); % Change (6,6) to any number you want

bestAIC = 99999;
bestm = [];
bestn= [];
for idx=1:npole
    for jdx=0:nzero
        if jdx>idx
            continue;
        end
        current_tf = tfest(z, idx, jdx);
        currentAIC = current_tf.report.Fit.AIC;
        fprintf('pole:%d, zero:%d, current AIC:%d, bestAIC: %d\n',idx,jdx,currentAIC,bestAIC);
        if currentAIC  < bestAIC
            bestAIC = currentAIC;
            bestm = idx;
            bestn = jdx;
        end
    end
end

besttf = tfest(z, bestm, bestn);
[n,d] = tfdata(besttf);
n = cellfun(@(x) {x.*(abs(x)>1e-7)}, n);
d = cellfun(@(x) {x.*(abs(x)>1e-7)}, d);
besttf = tf(n, d);

plot(z);
figure();
plot(0:z.Ts:z.Ts*(length(z.y)-1),z.y);
hold on;
[y2,t2]=lsim(besttf,z.u,0:z.Ts:z.Ts*(length(z.y)-1));
plot(t2,y2);
legend('original','systemID');
title('System ID result');

besttf

end