clc;
clear all;
close all;

%Zvolený prenos (volí sa taký, ktorý sa dá rozkmitať)
num = [1]; %volený čitateľ prenosu (numerator)
den = [5 14 7 1]; %volený menovateľ prenosu (denominator)
Gs = tf(num,den); %prenos systému
step(Gs,'r'); %prechodová charakteristika zvoleného systému
hold on;

%Výpočet kritických parametrov podľa metódy Zieglera a Nicholsa
[Gm,Pm,Wcg,Wcp] = margin(Gs); 
%Gm - Gain margin (bezpečnosť voči zosilneniu)
%Pm - Phase margin (fázová bezpečnosť)
%Wcg - The frequency where the gain margin is measured
%Wcp - The frequency where the phase margin is measured

if isinf(Gm)
    Gr = tf(0);
    disp('Sústavu sa nedá rozkmitať, Ziegler-Nicholsova metóda zlyhala');
    return
end

r0krit = Gm %kritické zosilnenie
Tkrit = 1/(Wcg/(2*pi)) %kritická perióda kmitov

%Vypočítané parametre klasického PID regulátoru podľa Zieglera a Nicholsa
r0 = 0.6*r0krit; %zosilnenie sústavy
Ti = 0.5*Tkrit; %integračná časová konštanta
Td = 0.12*Tkrit; %derivačná časová konštanta

%Prenos PID regulátoru
GrP = tf(r0);
GrI = tf(r0,[Ti 0]);
GrD = tf([r0*Td 0],1);
GrPID = GrP+GrI+GrD;

%Prechodová charakteristika regulačného obvodu 
title('Prechodová charakteristika systému');
xlabel('Čas t [s]');
ylabel('h(t)');
axis([0 80 -0.1 1.9]);
grid on;
regulacni_obvod=feedback(series(GrPID, Gs), 1);
step(regulacni_obvod,'b-');
legend('P.CH. systému bez regulátoru','P.CH. systému s PID regulátorom');