clear all; clc; close all;
%  Ajuste según Ziegler-Nichols
% tanque calefaccionado

pkg load control
s = tf('s');
t = 0:0.05:30;

% PROCESO
K1 = 1; tau1 = 5; K = 1; tau = 5; tita =  0.5;
p1 = K1/(tau1*s+1); p = K/(tau*s+1);
% uso aproximación de Padé de 2º orden
p2 = p*(1-tita/2*s+tita^2/12*s^2)/(1+tita/2*s+tita^2/12*s^2);

% TERMOCUPLA  TRANMITTER
Km = 0.32; taum = 0.1;
TT = Km/(taum*s+1);

% TRANSDUCTOR I/P
KIP = 1;

% VALVULA DE CONTROL
Kv = 1; tauv = 2;
CV = Kv/(tauv*s+1);

% CONTROLADOR PI
Kc = 1;
tauI = 100;
CON = Kc*(1+1/(tauI*s));

% bucle de control al aplicar escalón en set point
NUM = CON*KIP*CV*p2;
DEN = 1 + CON*KIP*CV*p2*TT;
rp = Km*NUM/DEN; 
rp= minreal(rp); %simplifica el numerador y el denominador
[num,den] = tfdata(rp, 'v');

y = step(rp,t);
plot(t,y)
xlabel('t')
ylabel('y')
title('Respuesta a escalón en set point - Ziegler-Nichols')
disp('Kc = ');disp(Kc)
disp('¿Se llegó a oscilación continua?')
h = input('¿si/no?  ','string');
while h == 'no'
  Kc = input('Kc = ');
    CON = Kc*(1+1/(tauI*s));
    NUM = CON*KIP*CV*p2;
    DEN = 1 + CON*KIP*CV*p2*TT;
    rp = Km*NUM/DEN; rp= minreal(rp); 
    [num,den] = tfdata(rp, 'v');

    y = step(rp,t);
    plot(t,y)
    xlabel('t');ylabel('y');
    title('Respuesta a escalón en set point - Ziegler-Nichols')
    disp('Kc = ');disp(Kc)
    h = input('¿si/no?  ','string');
end
if h == 'si' 
   Kc = input('Kcu = ');
    CON = Kc*(1+1/(tauI*s));
    NUM = CON*KIP*CV*p2;
    DEN = 1 + CON*KIP*CV*p2*TT;
    rp = Km*NUM/DEN; rp= minreal(rp); 
    [num,den] = tfdata(rp, 'v');

    y = step(rp,t);
    plot(t,y)
    xlabel('t');ylabel('y');
    title('Respuesta a escalón en set point - Ziegler-Nichols')
end

Pu = input('Pu = ')
Kc = 0.45*Kc; tauI = Pu/1.2;
CON = Kc*(1+1/(tauI*s));
NUM = CON*KIP*CV*p2;
DEN = 1 + CON*KIP*CV*p2*TT;
rp = Km*NUM/DEN; rp= minreal(rp); 

y = step(rp,t);
plot(t,y)
xlabel('t')
ylabel('y')
title('Respuesta a escalón en set point - Ziegler-Nichols')

disp(Kc); disp(tauI);