Choosing best system identification model using AIC

findbest.m

clear; close all;

s=tf('s');
G=1/(s^2+s+1);
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); %uを転置してカラムベクトルに直す
plot(z);

Options = tfestOptions;      
Options.Display = 'off';      
Options.WeightingFilter = [];

bestAIC = 99999;
m = 6;
n = 6;
bestm = [];
bestn= [];
for idx=1:m
    for jdx=0:n        
        if jdx>idx 
            continue;
        end
        tf = tfest(z, idx, jdx, Options);
        currentAIC = 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, Options)
figure();title('System ID result');
[y,t]=lsim(G,u,t);
plot(t,y);
hold on;
[y2,t2]=lsim(besttf,u,t);
plot(t2,y2);
legend('original','systemID');