%% Control por prealimentación --feedforward-- en proceso con retardo en perturbación.
% © 2020, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/ffret.html
% Este código funcionó con Matlab R2020b
% Objetivo: Comparar opciones de prealimentación en tiempo (step) y frecuencia (Bode) en un caso de estudio que incorpora retardo en una de las FdTs.
%
% Tabla de Contenidos
% Modelo
% Cálculo de prealimentación
%     Opción 1: estática
%     Opción 2: dinámica (no vale)
%     Opción 2: dinámica (con polo rápido para realizabilidad)
%     Opción 3: anticipar el control
%     Comparación de las diversas opciones en tiempo y frecuencia (sólo prealim.)
%% Modelo
% y=Gu+G_{pert}d

s=tf('s');

G=2/(s+1)^2;

G_pert_sinretardo=1/(2*s+1); retpert=1.2;
G_pert=exp(-retpert*s)*G_pert_sinretardo

step(G,G_pert), legend , grid on

%% Cálculo de prealimentación
%% Opción 1: estática

FF1=-tf(dcgain(G)\dcgain(G_pert))

%% Opción 2: dinámica (no vale)

FF2=-inv(G)*G_pert %NO realizable

%% Opción 2: dinámica (con polo rápido para realizabilidad)

FF2=FF2/(0.4*s+1); %jugar con este polo, para ver efecto sobre saturación/alta frecuencia en U. Diseño final, 0.4
zpk(FF2)

% y=GU_{FF}+G_{pert}d=(G*FF+G_{pert})d

step(G*FF1+G_pert,G*FF2+G_pert), grid on, legend("estática","dinámica")
step(FF1,FF2), grid on, legend("estática","dinámica")
bodemag(FF1,FF2), grid on, legend("estática","dinámica")

%% Opción 3: anticipar el control
% Como un retardo se puede aproximar por e^{-ds}=\frac{1}{e^{ds}}\approx \frac{1}{1+ds}, podemos hacer esa aproximación con las "dos terceras partes" del retardo:

retpert*0.333 %coincide con la cte. de tiempo del polo adicional de opc.2
FF3=-inv(G)*G_pert_sinretardo/(retpert*0.333*s+1)*exp(-retpert*0.667*s);
zpk(FF3)

%
%% Comparación de las diversas opciones en tiempo y frecuencia (sólo prealim.)
% Esta sería la acción ante escalón en perturbación

step(FF1,FF2,FF3), legend, grid on, yline(0,'-.k'), %ylim([-2 0.2])
title("Acción de control para cancelar un escalón de perturbación")
legend("Estático","Dinámico polo rápido","Dinámico v2","location","best")
bodemag(FF1,FF2,FF3), grid on
legend("Estático","Dinámico polo rápido","Dinámico  v2","location","best")
title("Acción de control (amplitud) para cancelar perturb. senoidal")

% Y el error residual sería:

yff1=G*FF1+G_pert;
yff2=G*FF2+G_pert;
yff3=G*FF3+G_pert;
step(yff1,yff2,yff3,G_pert), grid on, legend("Estático","Dinámico polo rápido","Dinámico  v2","sin prealim.")
title("Salida compensada con prealimentación ante pert. escalón")
bodemag(yff1,yff2,yff3,G_pert,logspace(-2,2,400)),ylim([-40 4]), grid on 
title("Bode de perturbación a salida")
legend("Estático","Dinámico polo rápido ","Dinámico v2","sin prealim.","Location","best")

%
%%