clc; clear
close all

s=tf('s'); 
t=0:0.01:1.2;

% Planta
Gp=4e5/(s*(s+20)*(s+200));

syslcc_nc=feedback(Gp,1);
figure; step(syslcc_nc);
axis([0 1.2 0 1.8])
% Access with a Structure field access (.)
spec=stepinfo(syslcc_nc);
%         RiseTime: 0.1833
%     SettlingTime: 2.8384
%      SettlingMin: 0.7211
%      SettlingMax: 1.5140
%        Overshoot: 51.3964
%       Undershoot: 0
%             Peak: 1.5140
%         PeakTime: 0.4854

tr_sc=spec.RiseTime;        % tiempo de subida sistema lazo cerrado sin compensación
ts_sc=spec.SettlingTime;    % tiempo de asentamiento sistema lazo cerrado sin compensación
tp_sc=spec.PeakTime;        % tiempo de pico sistema lazo cerrado sin compensación
Mp_sc=spec.Overshoot;       % sobrepaso sistema lazo cerrado sin compensación

% Compensador de adelanto diseñado en tiempo continuo
Kc=0.55;
Gc=Kc*(s + 15)/(s + 40);

Gla=Gc*Gp;
gan_la_sc=dcgain(Gla);
Glc=feedback(Gla,1);

spec=stepinfo(Glc);
tr_cc=spec.RiseTime;        % tiempo de subida sistema lazo cerrado con compensación

figure; step(syslcc_nc,Glc)
title('Desempeño del Sistema en Tiempo Continuo')
legend('Sin Compensar','Compensado');

%% Determinación del periodo de muestreo en base a las especificaciones 
% del sistema de control en lazo cerrado sin compensación
% Nm = tr/T
Nm=8;              % Número de muestras
T=tr_cc/Nm;        % Periodo de muestreo en base al tiempo de subida

% T=0.08;
z=tf('z',T);
Gpd=c2d(Gp,T,'zoh');
% figure; pzmap(Gpd)

%  Transfer function:
%         0.0009953*(z + 3.3503)*(z + 0.2396)
% Gpd = ---------------------------------------
%          (z - 1)*(z - 0.8973)*(z - 0.7177)
     
% syslcd_nc=feedback(Gpd,1);
% figure; step(syslcc_nc,syslcd_nc);
% title('Desempeño del Sistema en Tiempo Continuo y Discreto')
% legend('Sin Compensar Continuo','Sin Compensar Discreto');

% close all

%% Discretizacion del controlador en tiempo continuo usando aproximación de Tustin
% para verificar el efecto de la compensación polo-cero
Gcd=c2d(Gc,T,'tustin');
epsilon=0.0;
Cz=2*(z-epsilon)/(z + 1 - 2*epsilon);
Dz=Gcd*Cz;

Glad_cpz=minreal(Dz*Gpd);
gan_la_sdcpz=dcgain(Glad_cpz);

% figure; pzmap(Glad);
% Verificar estabilidad
figure; margin(Glad_cpz)
legend('Estabilidad del sistema en lazo cerrado con compensacion del ZOH');

Glcd_cpz=feedback(Glad_cpz,1);
figure; step(Glc,Glcd_cpz);
title('Desempeño del Sistema en Tiempo Continuo y Discreto Compensados')
legend('Sistema compensado tiempo continuo','Sistema compensado tiempo discreto CC Polo-Cero');

% Verificar estabilidad y desempeño del sistema discreto en lazo cerrado
% sin compensación del ZOH
Glad_scpz=Gcd*Gpd;
gan_la_sdspz=dcgain(Glad_scpz);

Glcd_scpz=feedback(Glad_scpz,1);
figure; step(Glc,Glcd_scpz);
title('Desempeño del Sistema en Tiempo Continuo y Discreto Compensados')
legend('Sistema compensado tiempo continuo','Sistema compensado tiempo discreto SC Polo-Cero');

% Verificar estabilidad
figure; margin(Glad_scpz); grid
legend('Estabilidad del sistema en lazo cerrado sin compensacion del ZOH');

figure; bode(Glad_cpz,Glad_scpz); grid
title('Comparacion de la estabilidad del sistema con y sin compensacion del ZOH');
legend('Sistema con compensacion del ZOH','Sistema sin compensacion del ZOH');

% Verificar desempeño en frecuencia del sistema en tiempo discreto compensado
% con y sin compensación del ZOH
figure; bode(Glcd_cpz,Glcd_scpz); grid
title('Desempeño en frecuencia del sistema discreto con y sin compensación del ZOH')
legend('Sistema con compensacion del ZOH','Sistema sin compensacion del ZOH');

% close all

%% Discretizacion del controlador en tiempo continuo usando aproximación de ZOH
% para verificar el efecto de la compensación polo-cero
% Gcd=c2d(Gc,T,'zoh');
% Cz=2*z/(z + 1);
% Dz=Gcd*Cz;
% 
% Glad_cpz=Dz*Gpd;
% gan_la_sdcpz=dcgain(Glad_cpz);
% 
% % figure; pzmap(Glad);
% % Verificar estabilidad
% figure; margin(Glad_cpz)
% legend('Estabilidad del sistema en lazo cerrado con compensacion del ZOH');
% 
% Glcd_cpz=feedback(Glad_cpz,1);
% figure; step(Glc,Glcd_cpz);
% title('Desempeño del Sistema en Tiempo Continuo y Discreto Compensados')
% legend('Sistema compensado tiempo continuo','Sistema compensado tiempo discreto CC Polo-Cero');
% 
% % Verificar estabilidad y desempeño del sistema discreto en lazo cerrado
% % sin compensación del ZOH
% Glad_scpz=Gcd*Gpd;
% gan_la_sdspz=dcgain(Glad_scpz);
% 
% Glcd_scpz=feedback(Glad_scpz,1);
% figure; step(Glc,Glcd_scpz);
% title('Desempeño del Sistema en Tiempo Continuo y Discreto Compensados')
% legend('Sistema compensado tiempo continuo','Sistema compensado tiempo discreto SC Polo-Cero');
% 
% % Verificar estabilidad
% figure; margin(Glad_scpz)
% legend('Estabilidad del sistema en lazo cerrado sin compensacion del ZOH');
% 
% figure; bode(Glad_cpz,Glad_scpz)
% title('Comparacion de la estabilidad del sistema con y sin compensacion del ZOH');
% legend('Sistema con compensacion del ZOH','Sistema sin compensacion del ZOH');
% 
% % Verificar desempeño en frecuencia del sistema en tiempo discreto compensado
% % con y sin compensación del ZOH
% figure; bode(Glcd_cpz,Glcd_scpz);
% title('Desempeño en frecuencia del sistema discreto con y sin compensación del ZOH')
% legend('Sistema con compensacion del ZOH','Sistema sin compensacion del ZOH');