%% Bucles de control ante referencia escalón, rampa, ...: ERROR estacionario; teoría + ejemplos Matlab
% © 2023, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Este código funcionó sin errores en Matlab R2022b (Linux)
% Objetivos: Comprender los conceptos de error de posición y velocidad, y desarrollar ejemplos Matlab de seguimiento de referencias escalón, rampa, parábola con distintos procesos y controladores. Las conclusiones NO son válidas ante perturbaciones a la entrada de integradores en el proceso.
% Presentaciones en vídeo:
% http://personales.upv.es/asala/YT/V/errposv1.html (teoría base)
% http://personales.upv.es/asala/YT/V/errposv2.html (error posición y velocidad)
% http://personales.upv.es/asala/YT/V/errposv3.html (entradas polinomiales)
% http://personales.upv.es/asala/YT/V/errposv4.html (ejemplos)
% Tabla de Contenidos
% Motivación
% Conceptos Teóricos
%         Error final ante referencia escalón (equilibrio)
%         Error final ante referencia rampa (rég. estacionario)
%     Errores de posición y velocidad según "tipo" de L(s)=G(s)K(s)
%         Tipo 0
%         Tipo 1
%         Tipo 2 o superior
%     Resultado teórico principal que engloba todo: error estacionario ante referencia polinomial
%         Resumen/Conclusiones
% Ejemplo 1: control P de sistema tipo 0 ante escalón y rampa de referencia
%         Referencia escalón
%         Referencia rampa
% Ejemplo 2: control PI de sistema tipo 0 ante escalón y rampa de referencia
%         Referencia escalón
%         Referencia rampa
% Ejemplo 3: control PID de sistema tipo 1 ante escalón y rampa de referencia
%         Referencia escalón
%         Referencia rampa
%
%% Motivación
% Referencia escalón es un cambio "brusco" de punto de funcionamiento; se desea saber si un sistema controlado se posiciona ahí (o cerca) cuando el transitorio acabe.
% Referencia "rampa" es un cambio "gradual" de punto de funcionamiento; se desea saber si es seguida con exactitud (error tiende a cero cuando acaba el transitorio), con retraso (misma pendiente, pero error constante en régimen permanente), o si no puede ser seguida (pendientes diferentes, error creciendo sin límite).
% En general, son casos particulares de referencias "polinomiales"... por ejemplo, seguir una "trayectoria de velocidad constante", una "trayectoria de aceleración constante" o similar puede ser un problema que se deba resolver en sistemas mecánicos/robótica, etc. Realmente más allá de la "aceleración" no tiene ya demasiado interés desde un punto de vista "general".
%
%
%
%% Conceptos Teóricos
% Diagrama de bloques de bucle cerrado (realimentación unitaria)
%
%
% El error de un bucle de control estándar es E(s)=\frac{1}{1+G(s)K(s)}r(s)=S(s)\cdot r(s). A la FdT S(s) se le suele llamar "función de sensibilidad" de bucle cerrado.
% También S(s) es la FdT entre "perturbación a la salida \delta_y" y la salida y o e (cambiando signo). Pero NO lo es entre "perturbación a la entrada \delta_u" y e, por lo que los resultados aquí presentados NO son válidos ante perturbación a la entrada con integradores en el proceso.
%
% Si denominamos L(s)=G(s)K(s), y expresamos en forma numerador/denominador L(s)=\frac{N_L(s)}{D_L(s)}, entonces
% S(s)=\frac{1}{1+\frac{N_L(s)}{D_L(s)}}=\frac{D_L(s)}{D_L(s)+N_L(s)}.
% A la FdT L(s) se le suele llamar función de bucle abierto. Obviamente, L(s) puede ser inestable ( D_L tener raíces de parte real cero o positiva), pero S ser estable (en este caso, raíces de D_L+N_L parte real negativa). Como S no tiene polos en cero, existe S(0).

syms s
G=3/(s+1)/s; K=2;
L=G*K
S(s)=simplify(1/(1+G*K)) %FdT de Referencia a Error de bucle cerrado

%
% Nota importante para ampliación/especialidad: control 2GL
% Aunque en este material se desarrolla el concepto de "error ante cambios de referencia" porque es la forma usual de explicarlo en una primera asignatura de control, realmente la importancia práctica podría estar más en el error ante perturbaciones de salida, porque los errores ante referencia son conocidos, se pueden precalcular, y por tanto puede modificarse la entrada al regulador para eliminarlo; eso es lo que se denomina control con 2 grados de libertad, discutido en los vídeos https://youtu.be/eBaleuoZe6M y https://youtu.be/Aq8sjhFX4ys , entre otros.
%
%% Error final ante referencia escalón (equilibrio)
% El teorema del valor final dice que, ante escalón unitario en referencia, el valor final del error valdrá (si el bucle cerrado es estable):
% e_p:=e(\infty)=\lim_{s\to 0} sE(s)=\lim_{s\to 0} s\cdot \left(S(s)\cdot \frac{1}{s}\right)=S(0)
% o sea, el error de posición es igual a la ganancia estática S(0) de la función de transferencia S(s), si S(s) es estable, obviamente.

S(0) %error FINAL ante escalón
ilaplace(S(s)*1/s) %error ante escalón

%% Error final ante referencia rampa (rég. estacionario)
% El teorema del valor final dice que, ante escalón unitario en referencia, el valor final del error valdrá (si el bucle cerrado es estable):
% e_v:=e(\infty)=\lim_{s\to 0} sE(s)=\lim_{s\to 0} s\cdot S(s)\cdot \frac{1}{s^2}=\lim_{s\to 0} \left(S(s)\cdot \frac{1}{s}\right)
% o sea, el error de velocidad e_v NO es igual a la ganancia estática de la función de transferencia S(s), sino a la ganancia de "S(s)/s", si existe. Claro, intuitivamente, la respuesta ante rampa de S(s) es la respuesta ante escalón de S(s)/s, por lo que le aplica lo visto arriba.
% Para que exista el límite de S(s)/s forzosamente S(s)=s\cdot \hat S(s), esto es, el numerador de S(s) tiene que tener una raíz al menos en s=0.

S(s)/s
limit(S(s)/s,0) %error FINAL ante rampa
ilaplace(S(s)*1/s^2) %error ante rampa

% Se le llama "error de velocidad" al valor final del error ante referencia rampa... aunque NO tiene dimensiones de "velocidad" (derivada de salida controlada), sino que tiene dimensiones de la propia "salida controlada". O sea, la "velocidad" a la que se refiere el término es a la de la referencia. Como puede dar lugar a cierta confusión, insistiremos más en los ejemplos posteriores en la interpretación correcta del término.
%
%% Errores de posición y velocidad según "tipo" de L(s)=G(s)K(s)
% El tipo de un sistema es el número de "integradores" en su FdT (número de polos en el origen).
% O sea, si L(s)=\frac{N_L(s)}{s^m\cdot \hat D(s)} donde ni N_L ni \hat D tiene raíces en cero, L(s) sería de "tipo m".
% El valor de " m " es la suma de los integradores que haya en el proceso G(s) + los integradores del controlador K(s).
% Obviamente, L(s) es una función inestable en bucle abierto si su tipo es mayor o igual a 1; obviamente, eso no supone ningún problema si en bucle cerrado se puede estabilizar.
% En bucle cerrado, con D_L(s)=s^m\cdot \hat D(s), con \hat D sin raíces en el origen, entonces
% S(s)=\frac{D_L(s)}{D_L(s)+N_L(s)}=\frac{s^m\cdot \hat D(s)}{s^m\cdot \hat D(s)+N_L(s)}
% el NUMERADOR de S(s) tiene m raíces en cero; el denominador D_L+N_L deberá tener todas sus raíces estrictamente en el semiplano izquierdo (parte real negativa) del plano complejo para poder continuar con la discusión, si no, el error tenderá a infinito incluso ante entradas acotadas.
%% Tipo 0
% S(s)=\frac{1\cdot \hat D(s)}{1\cdot \hat D(s)+N_L(s)}
% Ante referencia escalón:
% e_p=S(0) \neq 0, porque suponemos que \hat D(0)\neq 0, no hay raíces en el origen de \hat D.
% Ante referencia rampa: el error se hace infinito, porque en régimen estacionario tendrá pendiente S(0)\neq 0.

G(s)=3/(s+1);K=3;
S(s)=simplify(1/(1+G*K))
S(0)
ilaplace(S(s)/s) %error ante escalón
ilaplace(S(s)/s^2) %error ante rampa

%% Tipo 1
% S(s)=\frac{s\cdot \hat D(s)}{s\cdot \hat D(s)+N_L(s)}
% Ante referencia escalón:
% e_p=S(0) = 0, porque D_L(0)=0\cdot \hat D(0)=0.
% Ante referencia rampa:
% e_v=\lim_{s\to 0} S(s)\cdot \frac{1}{s}=\lim_{s\to 0}\left(\frac{s\hat D(s)}{s\hat D(s)+N_L(s)}\cdot \frac{1}{s}\right)=\frac{\hat D(0)}{N_L(0)}\neq 0
% Si el bucle es estable, entonces el error NO se hace infinito, tiende a un valor final (no cero).
% Resumen: siguen sin error referencias escalón, referencia rampa con error finito.
% Por superposición, referencia "rampa + escalón" tendría error finito, obviamente, sólo dependiendo de la "pendiente".

G(s)=3/(s+1)/s; K=3;
S(s)=simplify(1/(1+G*K))
S(0)
ilaplace(S(s)/s)%error ante escalón
limit(S(s)/s,0) %error de velocidad
ilaplace(S(s)/s^2) %error ante rampa

%%
%% Tipo 2 o superior
% e_p=0
% e_v=\lim_{s\to 0} S(s)\cdot \frac{1}{s}=\lim_{s\to 0}\frac{s^m\hat D(s)}{s^m\hat D(s)+N_L(s)}\cdot \frac{1}{s}=\lim_{s\to 0}\frac{s^{m-1}\hat D(s)}{N_L(s)}=0, para todo m>1.
% Resumen: siguen sin error referencias escalón o rampa, o combinaciones r(t)=a+b\cdot t.

G(s)=3/(s+1)/s; K=2*(s+1)*(s+3)/s;
S(s)=simplify(1/(1+G*K))
S(0)
simplify(ilaplace(S(s)/s),400) %error ante escalón
simplify(ilaplace(S(s)/s^2)) %error ante rampa

%
%% Resultado teórico principal que engloba todo: error estacionario ante referencia polinomial
% Un bucle de control estable donde L(s) sea de tipo "m" puede seguir sin error estacionario referencias r(t) que tengan forma de POLINOMIO de grado m-1.
% Si la referencia es un polinomio de grado m, el bucle cerrado tiene un error finito constante.
% Demostración:
% En Laplace, la referencia r(t)=a_0+a_1t+a_2t^2+\dots+a_{m-1}t^{m-1} tendrá una transformada dada por:
% r(s)=\frac{a_0}{s}+\frac{a_1}{s^2}+\frac{a_2\cdot 2!}{s^3}+\dots +\frac{a_{m-1}\cdot (m-1)!}{s^{m}}=\\
% \frac{a_0s^{m-1}+a_1s^{m-2}+2a_2s^{m-3}+\dots (m-1)!a_{m-1}}{s^m}=\frac{Q(s)}{s^m}
% Entonces, en efecto, el error de bucle, su valor final, será:
% e(\infty)=\lim_{s\to 0}s\cdot \left(\underbrace{S(s)\cdot r(s)}_{E(s)}\right)=\lim_{s\to 0}s\cdot \left( S(s)\cdot \frac{Q(s)}{s^m}\right)=\lim_{s\to 0} S(s)\cdot \frac{Q(s)}{s^{m-1}}
% Sustituyendo S(s), tendríamos:
% e(\infty)=\lim_{s\to 0} \frac{s^m\hat D(s)}{s^m\hat D(s)+N_L(s)}\cdot \frac{Q(s)}{s^{m-1}}=\lim_{s\to 0} \frac{s\hat D(s)Q(s)}{s^m\hat D(s)+N_L(s)}=0
% Obviamente, si fuera un polinomio de grado "m" la referencia, entonces r(s)=\frac{Q(s)}{s^{m+1}}, con lo que:
% e(\infty)=\lim_{s\to 0} \frac{s^m\hat D(s)}{s^m\hat D(s)+N_L(s)}\cdot \frac{Q(s)}{s^{m}}=\lim_{s\to 0} \frac{\hat D(s)Q(s)}{s^m\hat D(s)+N_L(s)} =  \frac{\hat D(0)Q(0)}{N_L(0)}\neq 0
%% Resumen/Conclusiones
% \bulletAnte referencia polinomio de grado "m-1", los sistemas de "tipo" m las siguen sin error. Si la referencia es un polinomio de grado m, la siguen con error finito constante. Si la referencia es un polinomio de grado \geq m+1, el error tiende a infinito.
% \bulletLo mismo dicho "a la inversa": si un sistema tipo m tiene error cero en régimen permanente ante cierta referencia, es porque el grado polinomial de la referencia es \leq m-1; si tiene error finito, no cero, es porque el grado polinomial de la referencia es m.
%
%
%
%
%% Ejemplo 1: control P de sistema tipo 0 ante escalón y rampa de referencia
% Pasemos a ejemplos Matlab:

syms s t
G(s)=8/(s+1); %es de tipo 0 (sin polos en el origen)
K=1.25;

% Sin acción integral en el control, el proceso es tipo 0, el controlador tipo 0, la función L(s) es de tipo 0.

S(s)=1/(1+G(s)*K)
S(s)=simplify(S(s))

%% Referencia escalón

e_p=S(0) %error de posición

% Es el error ante referencia constante unitaria (si el bucle cerrado es estable, obviamente). En efecto

err_s_step=S(s)*1/s
err_t_step=ilaplace(err_s_step)

% Su valor con t\to\infty coincide, claro, con teorema de valor final:

limit(s*err_s_step,0) %Resulta en E(0), la FdT evaluada en s=0

% La gráfica de la salida es:

fplot([1-err_t_step,1],[0 1],LineWidth=2), grid on, legend("salida control P","referencia",Location="best")

%% Referencia rampa
% Ante entrada "rampa" no es capaz de seguirla con la misma pendiente:

err_s_ramp=S(s)*1/s^2
err_t_ramp=ilaplace(err_s_ramp) %el error se hace "infinito" con el tiempo...

% El error aumenta, con pendiente 1/11, ganancia de S(s).
% El teorema de valor final da un resultado que realmente es inválido porque el límite de err_t_ramp no es finito:

limit(s*err_s_ramp,0) 

% Representemos la salida ante rampa

T1=G(s)*K/(1+G(s)*K);
T(s)=simplify(T1) %ref. a salida
simplify(T(s)==1-S(s)) %por definición e=r-y, entonces S(s)·r(s)=1·r(s)-T(s)·r(s)
y_t_ramp=ilaplace(T(s)*1/s^2)
fplot([y_t_ramp;t],[0 2],LineWidth=2), grid on, legend("salida control P", "referencia")

%
%% Ejemplo 2: control PI de sistema tipo 0 ante escalón y rampa de referencia

K=1.375*(s+1)/s; %PI por cancelación para tener los mismos polos

% Con acción integral en el control, el proceso es tipo 0, el controlador tipo 1, la función L(s) es de tipo 1.

S(s)=1/(1+G(s)*K)
S(s)=simplify(S(s))

%% Referencia escalón

e_p=S(0) %error de posición

% Es el error ante referencia constante unitaria (si el bucle cerrado es estable, obviamente). En efecto

err_s_step=S(s)*1/s
err_t_step=ilaplace(err_s_step)

% Su valor con t\to\infty coincide, claro, con teorema de valor final:

limit(s*err_s_step,0) %Resulta en E(0), la FdT evaluada en s=0
fplot([1-err_t_step,1],[0 1],LineWidth=2), grid on, legend("salida control PI","referencia")

%% Referencia rampa
% Ante entrada "rampa" sí es capaz de seguirla con la misma pendiente, y el error no se hace "infinito", aunque no se hace "cero":

err_s_ramp=S(s)*1/s^2
err_t_ramp=ilaplace(err_s_ramp) 

% Su valor con t\to\infty coincide, claro, con teorema de valor final:

e_v=limit(s*err_s_ramp,0) 

% *Nota: recordad que aquí el concepto de "ganancia estática" de S(s) no sirve para calcular este error e_v "error de velocidad", sino la ganancia de S(s)/s:

tmp_e_v(s)=S(s)/s
tmp_e_v(0) %este sí es el error de velocidad

% Dibujemos la salida, para comprobar:

T(s)=simplify(G(s)*K/(1+G(s)*K)) %ref. a salida
y_t_ramp=ilaplace(T(s)*1/s^2)
fplot([y_t_ramp;t],[0 2],LineWidth=2), grid on, legend("salida control PI", "referencia",Location="best")

% Tienen la misma pendiente, el error es "constante" al final del transitorio, porque T(s) tiene ganancia estática 1, y S(s) tiene ganancia estática cero: recordemos que la respuesta estacionaria ante rampa tiene de pendiente "la ganancia por la pendiente de la entrada", aunque estará desviada por una cierta constante...
% Cuando en bucles de control se habla de "error de velocidad", es importante darse cuenta que la velocidad es la correcta, exactamente igual a la de la rampa de referencia, o sea el "error de velocidad" sigue teniendo dimensiones de "posición".
%
%% Ejemplo 3: control PID de sistema tipo 1 ante escalón y rampa de referencia

G(s)=2/(s+1)/s;
K=4*(s+2)*(s+1)/s; %PID por cancelación parcial.

% Con acción integral en el control, el proceso es tipo 1, el controlador tipo 1, la función L(s) es de tipo 2.

S(s)=1/(1+G(s)*K)
S(s)=simplify(S(s)) % de ref. a error
T(s)=simplify(1-S(s)) %de ref. a salida

%% Referencia escalón

e_p=S(0) %error de posición

% Es el error ante referencia constante unitaria (si el bucle cerrado es estable, obviamente). En efecto

err_s_step=S(s)*1/s
err_t_step=ilaplace(err_s_step);

% Coincide, claro, con teorema de valor final:

limit(s*err_s_step,0) %Resulta en E(0), la FdT evaluada en s=0
fplot([1-err_t_step,1],[0 2.5],LineWidth=2), grid on, legend("salida control PI","referencia",Location="best")

%% Referencia rampa
% Ante entrada "rampa" sí es capaz de seguirla con error de velocidad cero:

err_s_ramp=S(s)*1/s^2
err_t_ramp=ilaplace(err_s_ramp) 

% Su valor con t\to\infty coincide, claro, con teorema de valor final:

e_v=limit(s*err_s_ramp,0) 

% Dibujemos la salida, para comprobar:

y_t_ramp=ilaplace(T(s)*1/s^2)
fplot([y_t_ramp;t],[0 2.5],LineWidth=2), grid on, legend("salida control PI", "referencia",Location="best")

% El sistema tipo 2 sigue sin error referencias polinomio grado 1, como la teoría predice.
% Seguirá con error finito referencias "parábola", aunque eso ya suele ser de menor interés práctico el subir grados:
% Dibujemos la salida, para comprobar:

y_t_parab=ilaplace(T(s)*2/s^3)
limit(y_t_parab-t^2,inf)
fplot([y_t_parab;t^2],[0 2],LineWidth=2), grid on, legend("salida control PI", "referencia",Location="best")

%