%% Análisis en frecuencia de bucle cerrado
% © 2019 Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en http://personales.upv.es/asala/YT/V/rfbcml.html
% Este código ha sido ejecutado en Matlab R2018b.
%% 1.- Proceso y Regulador

s=tf('s');
G=5/(s+1)^2/(0.01*s+1);

K1=.5*(1+.25*s/(0.05*s+1)+0.5/s)*1/(0.04*s+1); %K1 es tipo PID
zpk(K1)
K2=3.4*(1+.15*s/(0.01*s+1)); %K2 es tipo PD
zpk(K2)

[T1,S1,KS1]=cierrabucle(G,K1);
[T2,S2,KS2]=cierrabucle(G,K2);

%
%% 2.- Respuesta temporal en bucle cerrado

step(G/dcgain(G),T1,T2)
legend('control bucle abierto','bc1','bc2','Location','best'), grid on
step(tf(1/dcgain(G)),KS1,KS2,1.5), grid on
title('Acción de Control')
legend('control bucle abierto','bucle cerrado','bc2')
axis([0 1.5 -4 6])
%errores de posición
dcgain(T1),dcgain(S1),dcgain(KS1)
dcgain(T2),dcgain(S2),dcgain(KS2)

%
%% 3.- Respuesta en frecuencia en bucle cerrado
% Diagrama de SENSIBILIDAD S=1/(1+GK)

w=logspace(-1.25,2.25,150);
bodemag(tf(1),S1,S2,w), grid on
legend('bucle abierto','bc1','bc2','Location','best')
title('Sensibilidad referencia/pert. -> error')
%picos de sensibilidad
[val1,frec1]=norm(S1,inf)
val1dB=20*log10(val1)
[val2,frec2]=norm(S2,inf)
val2dB=20*log10(val2)
%simulemos en tiempo el peor caso de error del diseño 1:
t=0:0.01:9;
r=sin(frec1*t);
DosCosas=[T1;S1]; DosCosas.OutputName={'SALIDA','ERROR'};
lsim(DosCosas,r,t), grid on

% Diagrama de SENSIBILIDAD COMPLEMENTARIA T=1-S=GK/(1+GK)

w=logspace(-1,2,150);
bode(G/dcgain(G),T1,T2,w), grid on
legend('bucle abierto','bc1','bc2')
title('Sens. complem. referencia->salida')
%picos de sensibilidad complementaria
[val1,frec1]=norm(T1,inf)
val1dB=20*log10(val1)
[val2,frec2]=norm(T2,inf)
val2dB=20*log10(val2)
%simulemos en tiempo el peor caso de amplitud de salida en caso 2:
t=0:0.01:8;
r=sin(frec2*t);
T2.OutputName={'SALIDA control 2'};
lsim(T2,r,t), grid on

% Diagrama de SENSIBILIDAD DE CONTROL KS=K/(1+GK)

bodemag(tf(1/(G)),tf(1/dcgain(G)),KS1,KS2,w), grid on
legend('ideal u=inv(G)*r','bucl. abto','bc1','bc2','Location','best')
title('Sens. control referencia->acc.control')
%picos de sensibilidad de control
[val1,frec1]=norm(KS1,inf)
val1dB=20*log10(val1)
[val2,frec2]=norm(KS2,inf)
val2dB=20*log10(val2)
%simulemos en tiempo la respuesta ante ruido de medida de alta frecuencia para caso 2:
t=0:0.001:0.8;
r=0.05*sin(200*t);
DosCosas=[T2;KS2]; DosCosas.OutputName={'Salida y','CONTROL u'};
lsim(DosCosas,r,t), grid on 

%
%% Superposición de S, T, GK

bodemag(S1,T1,G*K1,w),grid on,
legend('S','T','GK')
title('diseño 1, superpuesto todo')
axis([0.1 100 -25 15])
bodemag(S2,T2,G*K2,w),grid on,
legend('S2','T2','GK2')
title('diseño 2, superpuesto todo')
axis([0.1 100 -26 20])

%

function [T,S,KS]=cierrabucle(G,K)
    T=feedback(G*K,1);%referencia a salida
    S=1-T; %referencia a error
    KS=minreal(K*S); %referencia a acción de control
end