%% Diseño de un IMC y un IMC-PID para un proceso de segundo orden + retardo [Simulación con SIMULINK]
% © 2019, Antonio Sala Piqueras. Universitat Politècnica de València.
% Este código ejecutó sin errores en Matlab R2019b
% Objetivos: Diseñar un IMC para modelo 2º orden+retardo, comprobar validez de especificaciones (no acelerar demasiado) y aproximarlo a un PID.
% Tabla de Contenidos
% 1. Modelado
% 2. Control IMC
%     Discusión sobre realizabilidad de Q (¿modelo de referencia orden 1?)
%     Discusión sobre validez de las especificaciones para un diseño IMC
% 3. Aproximación a PID (IMC-PID) en forma paralela
% Conclusiones
%%
%% 1. Modelado

s=tf('s');
tau1=4.2;tau2=2.8;
Modelo_sin_retardo=3/(tau1*s+1)/(tau2*s+1);
d=1;
Modelo=Modelo_sin_retardo*exp(-d*s);
zpk(Modelo)

%%
%% 2. Control IMC

Modelo_ref=1/(1.3*s+1)^2;
Q=inv(Modelo_sin_retardo)*Modelo_ref;
step(Modelo*Q,Modelo/dcgain(Modelo))
zpk(Q)

% Q sale realizable (2 polos y 2 ceros). Puede ser implementado explícitamente.
%% Discusión sobre realizabilidad de Q (¿modelo de referencia orden 1?)
%  Si se eligiera otra opción con modelo de referencia de orden 1:

Q_alt=inv(Modelo_sin_retardo)*1/(1.3*s+1);
zpk(Q_alt)

% El Q resultante tendría un cero más que polos.  Al transformar a PID (ver siguiente sección) daría un PID sin filtro de ruido. Ello requeriría algún tipo de filtro de alta frecuencia para hacerlo realizable (y limitar amplificación de ruido en la práctica).
%
%% Discusión sobre validez de las especificaciones para un diseño IMC
% Si acelerásemos más de 4 veces, grosso modo, la dinámica dominante, la respuesta ante perturbaciones de entrada no sería satisfactoria debido a la cancelación. Siguiendo el criterio de Skogestad como el polo dominante tiene constante de tiempo \tau_1=4.2, y nuestro modelo de referencia busca \tau_c=1.3, tenemos que 4\tau_c+\theta=6.2, con lo que esa constante de tiempo dominante \tau_1 es lo suficientemente rápida: las fórmulas del S-IMC harán que \min(\tau_1,4\tau_c+d)=\tau_1, por lo que el diseño de los PIDs de Skogestad S-IMC coincidirá con el IMC sin alteraciones.
% *Si \tau_1 hubiera sido mayor de 6.2, el proceso sería demasiado lento: no usaríamos el IMC literalmente para calcular un PID sino las modificaciones S-IMC (que, informalmente, diseñan la acción integral por lugar de las raíces).
%
%% 3. Aproximación a PID (IMC-PID) en forma paralela
% Construimos el regulador IMC implícito:

K_implicit_RetardoInterno=feedback(Q,Modelo,+1);

% Aproximamos retardo para tener un PID por serie de Taylor:
% e^{-ds}=1-ds+\frac{1}{2}d^2s^2-\frac{1}{6}d^3s^3+\dots

AproxRetardo=(1-d*s);
Modelo_SIMC=Modelo_sin_retardo*AproxRetardo;

% Construimos el regulador IMC implícito aproximado (debido al Taylor; sin aproximación tendría retardo interno):

K_implicit=minreal(feedback(Q,Modelo_SIMC,+1),1e-6); 
%minreal se para cancelar polos y ceros cercanos, tolerancia 1e-6

% El resultado es el regulador:

zpk(K_implicit) %es un PID con filtro de ruido de 2.13 rad/s de frec. de corte

% Si dispusiéramos de software que implementara controladores con funciones de transferencia arbitrarias, entonces ya habríamos terminado. Si sólo tenemos implementaciones de reguladores PID en determinadas parametrizaciones, debemos continuar.
% Descomponemos en fracciones simples para luego transformar en "parallel form PID" de Simulink

[r,p,h]=residue(K_implicit.num{:},K_implicit.den{:})

% Regulador expresado como:
% K_{implicit}=\frac{r_1}{s-p_1}+\frac{r_2}{s}+h=\frac{r_2}{s}+\frac{h\cdot s+(r_1-p_1h)}{s-p_1}
% Deseamos escribirlo como:
% K_{implicit}=K_d\cdot \frac{N\cdot s}{s+N}\ +\ \frac{K_i}{s}\ +\ K_p\ =\ \frac{K_i}{s}\ +\ \frac{(K_p+K_dN)\cdot s+K_pN}{s+N}

K_I=r(2)
N=-p(1)
syms Kp Kd real
sol=solve(Kp+Kd*N==h,Kp*N==r(1)-p(1)*h);
K_D=eval(sol.Kd)
K_P=eval(sol.Kp)
PID=K_D*N*s/(s+N)+K_I/s+K_P;
zpk(PID)

% Obviamente, no usaríamos diretamente la FdT que hemos llamado PID, sino que introduciríamos los parámetros K_D, K_P, K_I, N en el software que simulara los PIDs.
%
% Diseño alternativo con modelo de referencia de orden 1 (no realizable):

K_implicit_alt=minreal(feedback(Q_alt,Modelo_SIMC,+1),1e-6); 
%minreal se para cancelar polos y ceros cercanos

% El resultado es:

zpk(K_implicit_alt) %es un PID, sin filtro de ruido --no realizable sin él--

% Descomponiento en fracciones simples:

[r,p,h]=residue(K_implicit_alt.num{:},K_implicit_alt.den{:})

% Podemos obtener los coeficientes del PID

PID_alt = r/s + h(1)*s + h(2);
zpk(PID_alt)

% En una aplicación real, le añadiríamos, por ejemplo, un filtro de ruido 10 veces más rápido que el polo no dominante del proceso que incluíamos en el modelo:

polo2=1/tau2; N_alt=10*polo2
Filtro=N_alt/(s+N_alt);

% Podemos filtrar todos los componentes (prop., int., der.):

PID_alt_f1 = PID_alt*Filtro;
zpk(PID_alt_f1)

% O decidir sólo filtrar la acción derivada (si el sensor es lo suficientemente "limpio" para que el ruido en actuador debido a la acción proporcional sea aceptable), resultando en:

PID_alt_f2 = r/s + h(1)*s*Filtro + h(2);
zpk(PID_alt_f2)

% Dependiendo de la calidad del sensor y de los resultados experimentales observados, podría decidirse un ancho de banda mayor o menor del filtro de ruido.
%%
%% Conclusiones
% Hemos diseñado un IMC (y comprobado que no aceleramos "demasiado" la respuesta para necesitar las variaciones S-IMC), y aproximado mediante serie de Taylor el retardo para expresarlo como un PID.
% Para el proceso de 2o orden+retardo, un modelo deseado de 2o orden+retardo resulta directamente en un PID con filtro de ruido (lo hemos en forma paralela, pero podría hacerse de cualquier otra forma equivalente). Si el modelo de referencia hubiera sido de primer orden, el PID resultante tendría acción derivada "pura", sin filtrar, no realizable: el filtro debería añadirse a posteriori.