%% ME 4501
%% Compensacion Dinamica  
%% rhp, 2018,2021


%% En este script graficaremos el LGR de la funcion G(s)=1/s(s+1)
%% Imponemos un criterio z=0.5 que limita el sobrepaso (overshoot) 
%% Para cumplir con el criterio hay que compensar con un D(s) = (s+2)
%% Este cero deforma el LGR y se cumple el criterio para un rango mayor de K 
%% Sin embargo la funcion D(s) amplifica el ruido de alta frecuencia
%% AL agregar un polo de alta frecuencia, el ruido es cancelado

%% Para analizar estos efectos usamos graficos de los espectros en frecuencia


w = linspace(1e-3,1e2,1e6);                  % Eje de frecuencias (rd/s)

%% Funcion de transferencia G(s) = 1/s(s+1)
disp('Funcion de transferencia G(s) = 1/s(s+1)');
figure(1); 
G  = 1./( (j*w).^2 + (j*w) );                % en Fourier


subplot(2,1,1); 
plot((w), (abs(G))); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('|G(j\omega)|','fontsize',14); 
set(gca,'fontsize',12);title('G(j\omega) vs \omega');
axis([0 1 0 1e3]);

subplot(2,1,2); 
plot(log10(w), log10(abs(G))); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('log_{10}|G(j\omega)|','fontsize',14); 
set(gca,'fontsize',12); title('log_{10}(G(j\omega)) vs log_{10}(\omega)');

disp('Continuar');
pause;

%% Funcion de transferencia D(s) = (s+2)
%% El compensador D(s) = (s+2) aporta un cero que deforma el LGR

disp('Funcion de transferencia D(s) = (s+2)');
figure(2); 
D = (j*w) + 2; 
subplot(2,1,1); 
plot(w, abs(D)); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('|D(j\omega)|','fontsize',14); 
set(gca,'fontsize',12);title('D(j\omega) vs \omega');
%axis([0 1 0 1e3]);

subplot(2,1,2); 
plot(log10(w), log10(abs(D))); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('log_{10}|G(j\omega)|','fontsize',14); 
set(gca,'fontsize',12); title('log_{10}(G(j\omega)) vs log_{10}(\omega)');

disp('Continuar');
pause;

%% Funcion de transferencia D(s) = 1/(s+p) Agregando un polo Rapido
%% El polo cancelará el ruido de alta frecuencia al atenuar las altas frec.

disp('Funcion de transferencia D(s) = 1/(s+p) Agregando un polo Rapido');
figure(3); 
D = 1./(j*w + 15); 
subplot(2,1,1); 
plot(w, abs(D)); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('|D(j\omega)|','fontsize',14); 
set(gca,'fontsize',12);title('D(j\omega) = (s+2)/(s+p) vs \omega');
%axis([0 1 0 1e3]);

subplot(2,1,2); 
plot(log10(w), log10(abs(D))); grid; 
xlabel('\omega (rad/s)','fontsize',14); ylabel('log_{10}|G(j\omega)|','fontsize',14); 
set(gca,'fontsize',12); title('log_{10}(D(j\omega))) vs log_{10}(\omega)');


disp('Continuar');
pause;

%% En Octave Online graficamos LGR compensado
disp('graficamos LGR compensado');

num = [1]; den=[1 1 0]; G = tf(num,den);                 % G(s) inicial           en Laplace

num = [1 2]; den=[1 1 0]; G1 = tf(num,den);              % G1(s)=G(s)*(s+2)       en Laplace
num = [1 2]; den=conv([1 15],[1 1 0]); G2 = tf(num,den); % G2(s)=G(s)*(s+2)/(s+p) en Laplace

rlocus(G1);  % Grafica LGR de la Planta compensada por  D(s) = (s+2)
hold on; 
rlocus(G2);  % Grafica LGR de la Planta compensada por  D(s) = (s+2)/(s+15)
axis([-5 1 -3 3 ]);