
clc
clear all
close all

datos=load('DATOS_TEMPORALES_CALENTAMIENTO.txt');

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

figure(1)
title('Datos experimentales')
subplot(2,1,1)
plot(potencia)
title ('potencia')

subplot(2,1,2)
plot(temperatura)
title ('temperatura')


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 primer orden con retardo
% La ganancia es la misma que para la del primer orden.
% Queda por calcular el retardo y la nueva constante de tiempo.

cteTiempo=0.63*(temperatura(nMuestras)-temperatura(instanteEscalon));
while temperatura(instante)<(cteTiempo+temperatura(instanteEscalon))
    instante=instante+1;
end
disp('Parámetros del sistema:')
cteTiempo=instante-instanteEscalon
amplitudEscalon=(potencia(nMuestras)-potencia(instanteEscalon));
ganancia=(temperatura(nMuestras)-temperatura(instanteEscalon))/amplitudEscalon

% El número 1,6 corresponde al instante de muestreo en que la temperatura
% es 1,6ºC mayor que en el instante Escalón. Por prueba y error se ajusta
% el posible retardo
instante=instanteEscalon;
while temperatura(instante)<(temperatura(instanteEscalon)+1.6)
    instante=instante+1;
end

retardo=instante-instanteEscalon
cteTiempo=cteTiempo-retardo

modelo=tf(ganancia,[cteTiempo,1])
% la funsión setpDataOptions solo esta disponible para Matlab 2013b
opciones = stepDataOptions('StepAmplitude',amplitudEscalon);
[y,t] = step(modelo,opciones);

figure(2)
hold off
plot(t+retardo,y+temperatura(instanteEscalon),'linewidth',2)
hold on
plot([0:(nMuestras-instanteEscalon)],temperatura(instanteEscalon:nMuestras),'r','linewidth',2)
title('Primer orden con retardo')