%% Diseño de PID's por IMC: modificación para plantas lentas/integradores (Skogestad S-IMC)
% © 2019, 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/simcteo.html ,
% http://personales.upv.es/asala/YT/V/simcml.html .
% *Este código ejecutó sin errores en Matlab R2019b
%
% Objetivos: comprender como transformar un regulador IMC para un proceso sencillo (primer o segundo orden con retardo pequeño) en una estructura tipo PID cuando hay polos lentos o integradores... que, básicamente, se resume en que hay que dejar de utilizar el paradigma IMC.
% Tabla de Contenidos
% Revisión teórica IMC-PID para procesos sencillos
% Modificación Skogestad-IMC, teoría
%     Justificación
%     Resumen
% Ejemplo numérico IMC-PID versus S-IMC
% Apéndice: funciones auxiliares
%
% Modelo de la planta a utilizar en el ejemplo:

s=tf('s'); tau1=1; k=8; theta=0.05;
tildeG=k/(tau1*s+1); Gplus=exp(-theta*s);
G=tildeG*Gplus

%% Revisión teórica IMC-PID para procesos sencillos
% Consideremos el proceso: G=\frac{k}{(\tau_1s+1)(\tau_2s+1)}e^{-\theta s}, que podemos factorizar como  \tilde G=\frac{k}{(\tau_1s+1)(\tau_2s+1)}, G^+=e^{-\theta s}.
% Si selecionamos \tilde M=\frac{1}{\tau_cs+1}, el regulador IMC resultante sería:
% K=\frac{\tilde G^{-1}\tilde M}{1-\tilde M G^+}=\frac{(\tau_1s+1)(\tau_2s+1)}{k}\cdot  \frac{1}{\tau_cs+1-e^{-\theta s}}
%
% Aproximación del retardo: Para evitar estructuras tipo "predictor de Smith" y mantener la estructura PID, si \theta es suficientemente pequeño, a bajas frecuencias e^{-\theta s}\approx 1-\theta s, con lo que 1-e^{-\theta s}\approx 1-(1-\theta s)=\theta s, esto es:
% K\approx\frac{(\tau_1s+1)(\tau_2s+1)}{k}\cdot  \frac{1}{(\tau_c+\theta)s} =\frac{\tau_1}{k(\tau_c+\theta)}\cdot \frac{\tau_1s+1}{\tau_1s}\cdot (\tau_2s+1)
% que ya tiene una forma tipo PID serie K_c \cdot \frac{\tau_Is+1}{\tau_I s} \cdot \frac{\tau_Ds+1}{\tau_Fs+1}, con K_c=\frac{\tau_1}{k(\tau_c+\theta)}, \tau_I=\tau_1, \tau_D=\tau_2, habiendo añadido un filtro de ruido con \tau_F suficientemente pequeño para no alterar la dinámica dominante, pero suficientemente grande para evitar que el ruido de alta frecuencia se transmita al actuador.
% Inconvenientes: El IMC no funciona del todo bien cuando se trata de acelerar "mucho" al proceso... ante referencia GQ=M, pero ante perturbación de entrada (1-M)*G tiene los polos lentos de G, en general.
%% Modificación Skogestad-IMC, teoría
% S. Skogestad, en  https://doi.org/10.1016/S0959-1524(02)00062-8  , propone usar IMC sólo para acelerar hasta 4 veces, esto es usar \tau_I=\tau_1 (cancelando el polo (\tau_1s+1)) sólamente si las especificaciones de control verifican \tau_c \geq \tau_1/4.
% Bueno, realmente propone usar IMC hasta que \tau_1\leq 4(\tau_c+\theta), pero grosso modo dado que \theta debe ser pequeño, significa eso.
% En caso contrario, si \tau_1> 4(\tau_c+\theta) propone poner el cero de la acción integral en la cuarta parte de los polos deseados:
% K_{SIMC}=\frac{\tau_1}{k(\tau_c+\theta)}\cdot \frac{4(\tau_c+\theta)s+1}{4(\tau_c+\theta)s}\cdot (\tau_2s+1)
%
%% Justificación
% Se justifica cambiando la metodología de IMC a otra basada en "lugar de las raíces" (asignación de polos).
% Suponiendo \theta\approx 0 , si los polos de bucle cerrado deseados 1/\tau_c son mucho más rápidos que (\tau_1s+1), entonces el polo en 1/\tau_1 puede aproximarse a un integrador... en efecto, a frecuencias \omega mayores de 1/\tau_1, la planta G(s)=\frac{k}{\tau_1s+1}\approx \frac{k}{\tau_1s}; estamos suponiendo, en efecto, que eso ocurre cerca del "ancho de banda de bucle cerrado objetivo 1/\tau_c", si \tau_c<<\tau_1.
% Ahora, un controlador PI dado por \frac{\tau_1}{k\tau_c}(1+\frac{1}{4\tau_cs})=\frac{\tau_1}{k}\frac{(4\tau_cs+1)}{4\tau_c^2s} dará una función de lazo L=GK=\frac{4\tau_cs+1}{4\tau_c^2s^2}, con una ecuación característica 1+L=0 que, después de quitar denominadores será 0=4\tau_c^2s^2+4\tau_cs+1=(2\tau_cs+1)^2.
% Se puede comprobar entonces que:
% GK/(1+GK)\cdot r=\frac{4\tau_cs+1}{(2\tau_cs+1)^2}\cdot r                G/(1+GK)\cdot d_u=\frac{\frac{4k\tau_c}{\tau_1}s}{(2\tau_cs+1)^2}\cdot d_u
% ninguna función de bucle cerrado conserva el polo "lento" porque no se ha hecho cancelación.
% En esencia, la propuesta ante procesos con integradores es "desintonizar" la respuesta ante referencia, poniendo sus polos el doble de lentos de lo que M=1/(\tau_cs+1) diría, para a cambio acelerar la respuesta ante perturbaciones a la entrada del proceso.
% *El bucle cerrado real, sin la aproximación de G(s)=\frac{k}{\tau_1s+1}\approx \frac{k}{\tau_1s} sería: L=GK=\frac{\tau_1}{4\tau_c^2}\cdot \frac{(4\tau_cs+1)}{s(\tau_1s+1)}, que daría lugar a la ecuación característica 4\tau_c^2\tau_1\cdot s^2+4\tau_c(\tau_c+\tau_1)\cdot s+\tau_1=0, con dos polos entre \frac{-1}{\tau_c} y \frac{-1}{4\tau_c}.

tauc=tau1/15;L=tau1*(4*tauc*s+1)/4/tauc^2/s/(tau1*s+1);
rlocus(L), xlim([-1/tauc-1 0])
axis equal
hold on, plot(-1/tauc,0,'rx');
T=1/(1+L);
p1=plot(real(eig(T)),imag(eig(T)),'g*'); hold off
legend(p1,'polos b.c.')
text(-1/tau1-0.27,0.13,'-1/tau_1'), 
text(-1/tauc-0.27,-0.18,'-1/tau_c'),  text(-0.25/tauc-0.33,-0.26,'-1/(4tau_c)')

%
% Si \tau_1>>\tau_c, entonces tendria dos polos en \frac{-1}{2\tau_c}, como se ha visto antes.  En el otro extremo, si \tau_1=4\tau_c, quedaría el IMC/cancelación  L=GK=\frac{1}{\tau_cs} , T=GK/(1+GK)=\frac{1}{\tau_cs+1}.
%
%% Resumen
% Con lo que la regla S-IMC es escoger \tau_I=\min(\tau_1,4(\tau_c+\theta)) de modo que:
%     a) En control "poco agresivo" o proceso rápido \tau_1\leq 4(\tau_c+\theta) se sigue una estrategia IMC pura.
%     b) En control "agresivo" o proceso cuasi-integrador \tau_1\geq 4(\tau_c+\theta), se sigue estrategia de asignación de polos/evitación de cancelación. Para \tau_1\approx 4(\tau_c+\theta) el bucle cerrado tendría un polo cerca de ((\tau_c+\theta)s+1)y otro cerca de (\tau_1s+1), casi cancelado en la respuesta ante referencia... para \tau_1>>4(\tau_c+\theta) el sistema tendría un par de polos cerca de (2(\tau_c+\theta)s+1).
%
%% Ejemplo numérico IMC-PID versus S-IMC
% Especificaciones de control:

tauc=tau1/4; %constante de tiempo deseada (retardo aparte, claro)
t_est_objetivo= tauc*4+theta %si gastara IMC con retardo interno (pred. Smith)

% Diseño IMC-PID:

KIMC=(tau1*s+1)/(k*(tauc+theta)*s);zpk(KIMC)

% Diseño S-IMC:

tauI=min(tau1,4*(tauc+theta));
K_SIMC=tau1/k/(tauc+theta) * (tauI*s+1)/(tauI*s);zpk(K_SIMC)

% Simulacion y comparativa:

[T1,GS1,KS1]=SimulaBucle(G,KIMC,[],"NoPlot");
[T2,GS2,KS2]=SimulaBucle(G,K_SIMC,[],"NoPlot");
polosdeseados_bc=1/tauc
eig(pade(T1,2))
eig(pade(T2,2))
figure
step(T1,T2), grid on, legend('PID-IMC','S-IMC'), title('Step Ref->y')
step(KS1,KS2), grid on, legend('PID-IMC','S-IMC'),  title('Step Ref->u')
step(GS1,GS2), grid on, legend('PID-IMC','S-IMC'), title('Step du->y')

%
%% Apéndice: funciones auxiliares

function [T,GS,KS]=SimulaBucle(G,K,Tiempo,sinplot)
    arguments
     G
     K
     Tiempo = []
     sinplot = ""
    end
    T=feedback(G*K,1);
    GS=feedback(G,K); 
    KS=feedback(K,G); 
    if(sinplot == "NoPlot")
        return
    end
    subplot(1,3,1)
    step(T,Tiempo), grid on, title('Referencia->salida')
    subplot(1,3,2)
    step(KS,Tiempo), grid on, title('Referencia->control')
    subplot(1,3,3)
    step(GS,Tiempo), grid on, title('Perturb_u->salida')
end