%% Modelado, identificación y control PID de sistema de primer orden + integrador (control de eje mediante par manipulado)
%%
%
%
% © 2025, Antonio Sala, Universitat Politècnica de València, todos los derechos reservados.
% Objetivo: detallar un caso de estudio de control de posición o velocidad de un eje, desde el modelado físico hasta los diferentes controladores que podrían ser necesarios según la variable controlada elegida y las especificaciones estáticas/dinámicas exigidas.
%
% Presentaciones en vídeo:
% https://personales.upv.es/asala/YT/V/pideje1.html [planteamiento, modelado]
% https://personales.upv.es/asala/YT/V/eje2step.html [análisis de respuesta temporal bucle abierto]
% https://personales.upv.es/asala/YT/V/eje3idba.html [identificación experimental, bucle abierto]
%
% Tabla de Contenidos
% Modelado, identificación y control PID de sistema de primer orden + integrador (control de eje mediante par manipulado)
%
%     Objetivo general, motivación
%     1.- Modelado
%         Ecuaciones de la Física
%         Modelo en el dominio de Laplace
%     2.- Análisis de la respuesta temporal
%         Velocidad ante escalón de par
%         Posición ante escalón de par
%     3.- Identificación experimental
%
%% Objetivo general, motivación
% Haremos un "ejercicio académico" para modelar, identificar experimentalmente y controlar velocidad o posición del sistema de la figura:
%
%
% El sistema tendrá dos entradas ( T_{mot}, T_{carga} ) y dos estados (posición y velocidad angular; 1 grado de libertad mecánico).
% Para control, consideraremos el par producido por el motor T_{mot} como entrada de tipo variable "manipulada" (acción de control), y el par T_{carga} como entrada de tipo "perturbación".
% *El análisis estará incompleto (por simplicidad e ignorancia del que escribe esto), en muchísimos aspectos:
%
% A) Cosas pendientes desde un punto de vista de "ingeniería de control":
% \bullet aquí visualizaremos únicamente la "variable de salida controlada", pero se necesitaría visualizar más cosas, como...
% \bullet la gráfica de la acción de control (para comprobar si se satura o no), y
% \bullet posiblemente el diagrama de Bode desde "ruido de medida" hasta "acción de control", para comprobar si hay una excesiva amplificación del ruido de alta frecuencia. Casi seguro que necesitamos "filtrar" algo.
% \bullet Compensación de la fricción mecánica: el rozamiento lineal es simplista... la fricción seca (Coulomb) si es elevada, puede causar problemas en el control de "posición". Igual hay que hacer T_{mot}="resultado_del_PID"+ \gamma\cdot \mathrm{sign}(\omega)... sobre todo en control de posición (en control de "velocidad" podría considerarse una "perturbación escalón").
% \bullet ¿La "perturbación de carga" es realmente un "par" o es la "conexión mecánica de otro sistema con más inercia" que va a hacer todo más lento o más oscilatorio y va a cambiar las FdTs del modelo?
%
% B) Aparte, necesitamos discutir:
% \bullet Tecnología mecánica: diámetro de ejes, material de los mismos, uniones entre elementos, fatiga, vibraciones, engranajes y reductoras...
% \bullet "tecnología de implementación del control"  (¿analógica con amplificadores operacionales vs. digital con software vs. puramente mecánica con muelles y amortiguadores en caso de "posición"?)
%  \bullet "tecnología de actuación" (¿cómo proporcionar el par necesario? Motor eléctrico de CC, de alterna, inversor, brushless, ... Control "en cascada" --electrónica interna de dicho motor--).
% \bullet "tecnología de sensorización": cómo medir posición o/y velocidad de forma económica y fiable.
% \bullet Esto es un "subsistema" que está moviendo "otras cosas"... ¿Para qué lo queremos? ¿Por qué fijamos unas u otras especificaciones?
% \bullet ¿Hay que gestionar "comunicaciones" con quien le da la referencia? ¿Interfaz de usuario? ¿Registro de datos para diagnóstico en línea o post-mortem en caso de fallos?
% \bullet Moviendo cosas igual necesitamos setas de emergencia, fotocélulas que activen barrera de seguridad, etc. También cuadro eléctrico, magnetotérmicos, diferencial...
% \bullet Análisis económico de todo esto en la fase de diseño/proyecto: viabilidad / tecnologías alternativas (motor paso a paso, posicionamiento neumático, robot...).
%
% C) Ante toda esta "complicación"... existen soluciones "llave en mano" de control de ejes, "plug and play", de modo que a nivel de usuario o ingeniero que no piensa en I+D, todo lo de arriba es "transparente" a nosotros, nos lo han resuelto los que venden el equipo u otra ingeniería.
%
%
%% 1.- Modelado
%% Ecuaciones de la Física
% Un modelo como el de la figura tendrá de ecuaciones de estado
% \frac{d\omega}{dt}=\frac{1}{J}\left(-b\omega+T_{mot}-T_{carga}\right)
% \frac{d\theta}{dt}=\omega
% Nota 1: Si sólo buscamos control de 'velocidad' y disponemos de un sensor de dicha velocidad, entonces no sería necesaria la segunda ecuación. Pero para control de 'posición' necesitaremos ambas ecuaciones de estado, obviamente.
% Nota 2: tecnológicamente, necesitaremos dar el par motor con alguna tecnología como, por ejemplo, un motor eléctrico configurado en "control de par" ("control esclavo" de corriente en "cascada", se supone que mucho más rápido que la dinámica que deseemos del sistema mecánico bajo control).
%% Modelo en el dominio de Laplace
% A partir de las ecuaciones transformadas (c.i. nulas):
% s\cdot \omega(s)=\frac{1}{J}\left(-b\omega(s)+T_{mot}(s)-T_{carga}(s)\right)
% s\cdot \theta(s)=\omega(s)
% podremos despejar:
% \omega(s)=\frac{1}{Js+b}\cdot (T_{mot}(s)-T_{carga}(s))
% de modo que en control de velocidad, el proceso tendría FdT: \frac{1}{Js+b}.
%
%
% En control de posición, tendríamos:
% \theta(s)=\frac{1}{s}\omega(s)=\frac{1}{s\cdot(Js+b)}\cdot (T_{mot}(s)-T_{carga}(s))
% con lo que la FdT sería \frac{1}{s(Js+b)}.

syms J b real
syms s %Laplace variable
syms t real %time
Gvel(s)=1/(J*s+b)
Gpos(s)=Gvel/s;
Jnumerico=0.1; bnumerico=0.11; %Para ejemplos concretos

%% 2.- Análisis de la respuesta temporal
%% Velocidad ante escalón de par
% Podría ser conveniente expresar la FdT de varias formas:
% G_{vel}(s)=\frac{1}{Js+b}=\frac{1/J}{s+b/J}=\frac{1/b}{(J/b)\cdot s+1}
% El "polo de bucle abierto" es polo=-b/J, valor absoluto de parte real \sigma=b/J.
% En la forma K/(\tau s +1), la "constante de tiempo" es la inversa del polo \tau=J/b, y K=1/b.
% \bullet La respuesta libre de la velocidad tendrá términos e^{-\frac{b}{J}\cdot t}, esto es e^{polo\cdot t}, equivalementemente e^{-\sigma t} o bien e^{-\frac{1}{\tau}t, según nos convenga.
% \bullet La ganancia estática es G_{vel}(0)=\frac{1}{b}, de modo que \omega_{eq}=\frac{1}{b}\cdot (T_{mot,eq}-T_{carga,eq}), si existen los valores finales en los pares.
% \bullet La respuesta en velocidad ante escalón de par será:

omega_esc(s)=Gvel*1/s
partf_omega=partfrac(omega_esc(s),s,FactorMode="full")
omega_esc_t=ilaplace(partf_omega)

% El valor final es 1/b, el tiempo de establecimiento 98% sale de igualar e^{-bt_{est}/J}=0.02, que resulta aproximadamente t_{est}=\frac{4}{b/J}=\frac{4}{\sigma}, o t_{est}=4\frac{J}{b}=4\tau.
% Si sustituimos valores numéricos:

omega_numerico=subs(omega_esc_t,{b,J},{bnumerico,Jnumerico});
taunumerico=Jnumerico/bnumerico
Knumerico=1/bnumerico
vpa(omega_numerico,3)
fplot(omega_numerico,[0 8],LineWidth=3), grid on
xline(taunumerico,'r',Label="tau")
yline(0.63*Knumerico,Label="63%",Color=[.5 .6 .7])
xline(3*taunumerico,'g',Label="3*tau (95%)")
xline(4*taunumerico,'g',Label="4*tau (98%)")
title("Velocidad (rad/s) ante escalón unitario de par")

%
%% Posición ante escalón de par
% Tenemos:
% G_{pos}(s)=\frac{1}{s\cdot (Js+b)}=\frac{1/J}{s\cdot (s+b/J)}=\frac{1/b}{s\cdot ((J/b)\cdot s+1)}
% Como tiene un polo en s=0, la respuesta libre tendrá e^{-b/J\cdot t} y también e^{0t}=1, términos constantes que no vuelven a cero. La respuesta forzada ante escalón tendrá términos en "rampa".
% Por tanto, el sistema es "inestable", cualquier mínima perturbación en T_{carga} lo hará girar sin parar, y no podrá ser utilizado en la práctica para posicionarse en algún sitio en concreto sin "control en bucle cerrado".
% Desde un punto de vista "cualitativo" no hay mucho más que decir.
% De todos modos, podemos calcular la respuesta analíticamente para comprobar todo lo dicho:

pos_esc(s)=Gpos(s)*1/s
partf_pos=partfrac(pos_esc(s),s,FactorMode="full")
pos_esc_t=ilaplace(pos_esc)

% El transitorio exponencial durará lo mismo que el de velocidad, y el régimen estacionario (cuando la exponencial termine) será una rampa:
% pos_{est}(t)=\frac{1}{b}\cdot (t-\frac{J}{b})=K\cdot (t-\tau).

pos_numerico=subs(pos_esc_t,{b,J},{bnumerico,Jnumerico});
vpa(pos_numerico,4)

% Observemos que 9.09\cdot (t-\tau)=9.09\cdot (t-0.9091)=9.09\cdot t-8.264.

fplot(pos_numerico,[0 6],LineWidth=3), grid on
hold on
fplot(1/bnumerico*(t-taunumerico),[0 6],':',LineWidth=2), hold off
yline(0),xline(taunumerico,'r',Label="tau")
xline(3*taunumerico,':',Label="3·tau",Color=[.5 .6 .7])
legend("Posición(t), radianes","Rég. estacionario",Location="best")
title("Posición ante escalón unitario de par")
fplot(atan2(sin(pos_numerico/2),cos(pos_numerico/2)),[0 4],LineWidth=2.5), grid on
title("Posición (rad) ante escalón 0.5 de par (sensor monovuelta)")

%%
%% 3.- Identificación experimental
% En K/(\tau s+1), o K/s/(\tau s+1), se ha marcado en las figuras dónde encontrar el valor de la constante de tiempo \tau si se dispusiera de datos de laboratorio experimentales.
% En cuanto a la ganancia estática (parámetro K):
% \bullet En caso de "velocidad ante escalón", la ganancia K sería "incremento de velocidad"/"amplitud escalón de entrada".
% \bullet En caso de ensayo "posición ante escalón", , la ganancia K sería "pendiente de posición"/"amplitud escalón de entrada". En el segundo caso, debería ser T_{carga}=0 para observar la gráfica mostrada arriba. Eso puede no ser tan fácil de realizar "en la práctica", porque las perturbaciones usualmente están "fuera de mi control": no puedo "manipularlas"... y si el sistema es inestable, las perturbaciones podrían hacer "girar al sistema de forma descontrolada"... identificar sistemas inestables y con perturbaciones es más difícil en la práctica.
%
% *Identificación de parámetros físicos: a partir de las gráficas se puede determinar K y \tau.
% Con ello, como K=1/b y \tau=J/b, podrían despejarse los valores de J y b identificados.
%
%
%%