clc
clear all
close all

% El calentamiento se inicia al aplicarse un escalon de potencia.

datos=load('DATOS_TEMPORALES_CALENTAMIENTO.txt');

tiempo=datos(:,1);
potencia=datos(:,2);
temperatura=datos(:,3);
nMuestras=size(temperatura,1);

figure(1)
subplot(2,1,1)
plot(potencia,'linewidth',2)
grid on
axis([-100 4000 0 60]);
title ('Datos experimentales Acción de Control.');
xlabel('Muestras. (Se toma una muestra cada segundo)');
ylabel('% Potencia. ( 0 % -> 0 mA)');
subplot(2,1,2)
plot(temperatura,'linewidth',2)
grid on;
axis([-100 4000 15 95]);
title ('Datos experimentales Temperatura Planta.');
xlabel('Muestras. (Se toma una muestra cada segundo)');
ylabel('Temperatura Planta');

instante=1;
while(potencia(instante)==potencia(instante+1))
    instante=instante+1;
end
disp('El escalon se aplica en el instante:')
instanteEscalon=instante;

% Identificación de un modelo de segundo orden con retardo
% La ganancia es la misma que para la del primer orden.
% Queda por calcular la nueva constante de tiempo.

amplitudEscalon=(potencia(nMuestras)-potencia(instanteEscalon));
ganancia=(temperatura(nMuestras)-temperatura(instanteEscalon))/amplitudEscalon;

instante=instanteEscalon;
t72=0.72*(temperatura(nMuestras)-temperatura(instanteEscalon));
while temperatura(instante)<t72
    instante=instante+1;
end

t72=instante-instanteEscalon;
L=-612; %Este valor está calculado con pruba y error hasta que le respuesta
   %del modelo se ha ajustado completamente a la respuesta real del proceso
Ti=0.8*(t72-L);
tz=0.4*t72+0.6*L;
yz=temperatura(instanteEscalon+round(tz));
Yz_YKA=yz/temperatura(nMuestras);
% con el valor que sale de la division se busca en la tabla de la figura
% 5.13 de las transparencias de Identificación de procesos y la z tiene
% este valor
z=0.66;
tau72=1.258;
Ti=(t72-L)/tau72;
T1=z*Ti
T2=Ti-T1

syms s
modelo=tf(ganancia,sym2poly((T1*s+1)*(T2*s+1)))
opciones = stepDataOptions('StepAmplitude',amplitudEscalon);
[y,t] = step(modelo,opciones);
 
figure(2)
hold off
plot(t,y+temperatura(instanteEscalon),'linewidth',2)
hold on
plot(0:(nMuestras-instanteEscalon),temperatura(instanteEscalon:nMuestras),'r','linewidth',2)
grid on;
axis([-100 4000 15 95]);
title('Respuesta Escalon entrada. Sistema 2º Orden Sobreamortiguado');
xlabel('Muestras. (Se toma una muestra cada segundo)');
ylabel('Temperatura Planta');