Digital Control Octave Scripts

System Identification Simulation

% Digital Control System Identification Example
pkg load control
clear all;
close all;
% Set up the system to fit to.
A=diag([-1,-2]) % Make sure the system is stable.
B=[1,1]'
C=[1,2]
%A=diag(-1)
%B=1
%C=1
sys = ss(A,B,C); % sys is the system to fit.
% Move to the sampled domain.
T=0.01; % T is the sampling interval.
N=10000; % N is the number of training points.
t=0:T:N*T;
sys_discrete = c2d(sys,T);
sys_discrete_tf = tf(sys_discrete)
if(!isstable(sys)) 
  error('System is not stable!\n') % Check for stability.
 end
if(!isctrb(sys)) % Check to make sure the system is contrololable,
  error('System is not controllable!\n')
end
if(!isobsv(sys)) % and observable.  If it is not you need to fix it!
  error('System is not observable!\n')
end
% Now set up the system identification.
u=randn(size(t));
[y,t]=lsim(sys,u,t);
% The order of the 0fitted system is n.
n=2
Y=[[y(n+1:N+1)]]; % Note: Octave's lowest index is 1, not 0.
for k=n:N
  PSI(k-n+1,:)=[y(k:-1:k-n+1)', u(k:-1:k-n+1)];
end
PSI;
THETA=(PSI'*PSI)^(-1)*PSI'*Y
sys_fit=tf(THETA(n+1:2*n),[1,-THETA(1:n)'],T) % tf(num,den,T)
[y_fit,t]=lsim(sys_fit,u,t);
err=(y-y_fit);
plot(t,err)
figure;
plot(t,y,t,y_fit)
sum_err = (y-y_fit)'*(y-y_fit)