clc;
clear all;
close all;
 
%Prenos systému 
num = [1]; %čitateľ prenosu (numerator) 
den = [5 14 7 1]; %menovateľ prenosu (denominator)
Gs = tf(num,den); %prenos systému
Ts = 0.01; %vzorkovacia perióda 
t = (0:Ts:150); %časový interval
 
%Simulácia prenosu
u = [0 ones(1,length(t)-1)];
y = lsim(Gs,u,t);
title('Časová odozva dynamického systému');
plot(t,y,'r'); %vykreslenie simulácie systému
hold on;
 
%---------------------------------------------------------------------
%Aproximácia systému pomocou MNŠ prenosom 2. rádu  
 
phi = [ ]; %matica nameraných údajov
 
%Výpočet neznámych parametrov pomocou metódy najmenších štvorcov
 
for i = 2:-1:1
   phi = [phi y(i:end-2+i-1)]; 
end
 
for j = 2:-1:1
   phi = [phi u(j:end-2+j-1)']; 
end
 
Mvyst = y(3:end);
theta = (phi'*phi)^-1*phi'*Mvyst; 
 
%Určenie prenosu systemu
GzIdent = tf([theta(3:4)'],[1 -theta(1:2)'],(Ts)) %diskrétný prenos
GsIdent = d2c(GzIdent) %prevod diskrétneho prenosu na spojitý
yIdent = lsim(GsIdent,u,t);
title('Časová odozva systému');
xlabel('Čas t [s]');
axis([0 50 -0.1 1.1]);
grid on;
plot(t,yIdent,'b'); %vykreslenie diskrétneho identifikovaného prenosu
legend('Prenos sústavy','Odhad prenosu sústavy modelom druhého rádu', 'Location', 'southoutside');