%% Respuesta ante ESCALÓN de procesos de SEGUNDO orden sin ceros, para análisis de respuesta o para identificación experimental
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentaciones en vídeo:
%  --- http://personales.upv.es/asala/YT/V/ord2moti.html
%  --- http://personales.upv.es/asala/YT/V/ord2step.html
% Este código funcionó sin errores con Matlab R2022b (Linux)
% Objetivo: comprender la teoría y la identificación experimental rápida en el caso de respuesta ante escalón de sistema lineal de segundo orden (sin ceros) tanto en el caso oscilatorio (subamortiguado) como en el NO oscilatorio (sobreamortiguado).
%
% Tabla de Contenidos
% Motivación
% Polos de un sistema de segundo orden en forma "coef. de amortiguamiento (xi), frecuencia natural (omega_n)"
% Respuesta escalon de sistemas de 2o orden oscilatorios sin ceros
%     Parámetros característicos de la respuesta
%         Valor final
%         Tiempo de establecimiento
%         Oscilaciones: frecuencia propia, tiempo de pico, valor máximo
% Ejemplo numérico
%         Otros parámetros de posible interés (tpo. retardo / subida)
%     Conclusión del ejemplo: descripción en términos "útiles" de la respuesta
%% Motivación
% En general, los temas de "análisis de la respuesta de sistemas sencillos" se motivan porque es de gran interés poder describir de forma "entendible" dicha respuesta a un usuario/cliente (una "sarta de exponenciales o/y fórmulas trigonométricas" asusta a mucha gente). A la inversa, observando la respuesta, es de interés tener reglas "rápidas" de identificación experimental.
% Los sistemas RLC, masa-muelle-amortiguador, etc. tienen una FdT de segundo orden que, puede ocurrir que:
% Tenga polos reales: transitorio no oscilatorio, por lo que sólo tendremos interés en el valor final (derivadas a cero, G(0)) y el "tiempo de establecimiento", igualando exponenciales a 0.05 o 0.02. Físicamente, ocurre con pérdidas elevadas en resistencia y amortiguadores. Este caso se suele denominar sobreamortiguado en los textos de dinámica de sistemas.
% Tenga polos complejos: transitorio oscilatorio. Tendremos interés en unas cuantas cosas más... ese es el objetivo de este material. Este caso se suele denominar subamortiguado en los textos.

G1=tf(36,[1 7.5 9])
p1=pole(G1) %reales
G=tf(36,[1 2.4 9])
pole(G) %complejos
step(G1,G), grid on, legend("transitorio no oscilatorio (sobreamort.)","oscilatorio (subamortiguado)")

%% Polos de un sistema de segundo orden en forma "coef. de amortiguamiento (xi), frecuencia natural (omega_n)"
% Un sistema de segundo orden lineal dado por la EDO:
% \frac{d^2y}{dt^2}+2\xi\omega_n\frac{dy}{dt}+\omega_n^2y=K\omega_n^2\cdot u
% con condiciones iniciales nulas, representado en el dominio de Laplace con  Y(s)=G(s)U(s), siendo
%  G(s)=\frac{K\omega_n^2}{s^2+2\xi\omega_n s +\omega_n^2}
%  su función de transferencia, con \omega_n>0, 0< \xi <1 [subamortiguado], tiene unos polos (raíces del denominador) dados por:
% \frac{1}{2}\left(-2\xi\omega_n\pm \sqrt{(2\xi\omega_n)^2-4\omega_n^2}\right)=
% \frac{1}{2}\left(-2\xi\omega_n\pm 2\omega_n\sqrt{\xi^2-1}\right)=\left(-\xi\pm j\cdot \sqrt{1-\xi^2}\right)\omega_n
% cuya parte real es -\xi\omega_n, y parte imaginaria \omega_p=\pm\omega_n \sqrt{1-\xi^2}. ESTABLE.
%
% *Nota 1: ¿por qué este G(s) y no G(s)=\frac{a}{s^2+bs+c}? Pues por razones "estéticas" y "mnemotécnicas": salen "bonitas" las fórmulas más abajo para memorizar, pero en esencia es obviamente lo mismo.
%
% *Nota 2: Comportamiento con otros valores de \xi:
% Si |\xi|>1 [sobreamortiguado] entonces {\xi^2-1}>0 y, por tanto \sqrt{\xi^2-1} es un número REAL. Salen "exponenciales no oscilatorias" que no son objetivo principal de aquí.
% Si \xi\leq 0 entonces  \sigma=-\xi\omega_n\geq 0, salen exponenciales "positivas" INESTABLES (oscilatorias si -1<\xi\leq 0, no oscilatorias si \xi\leq -1). Los procesos "inestables" son de poco interés (en principio): no pueden ponerse en marcha porque algo se rompe.
% Ejemplo:

G1
K1=3; omega_n=3; xi= 7.5/2/omega_n
p1=pole(G1)
log(0.02)./p1 %tiempo de establecimiento 98%, lo único que hay que decir en caso sobreamortiguado.
4./abs(p1) 
G
K1=3; omega_n=3; xi= 2.4/2/omega_n
pole(G)

%
%% Respuesta escalon de sistemas de 2o orden oscilatorios sin ceros
%
% Calculemos la respuesta temporal ante u(t)=A, o sea  U(s)=A\cdot \frac{1}{s} en dominio de Laplace, dada por la fórmula que obtendrá la Symbolic toolbox con:

syms A K xi omega_n real
assume(0<=xi<1) %par de raíces complejas
assume(omega_n>0) %si no, saldrá inestable...
syms s
Y= A/s   *   K*omega_n^2/(s^2+2*xi*omega_n*s+omega_n^2)
syms t real
Y(t)=simplify(ilaplace(Y))

% Lo "importante" es que son senos y cosenos modulados por exponencial que se suman al punto de equilibrio (valor final) AK.
%% Parámetros característicos de la respuesta
%% Valor final
% Es AK. También sale de la ecuación y_{eq}=G(0)u_{eq}, válida para cualquier FdT.
%% Tiempo de establecimiento
% El transitorio dura (igualando la exponencial a 0.05), criterio 95%

vpa(log(0.05)/-omega_n/xi,3)

% También, el transitorio dura (igualando la exponencial a 0.04321), criterio 95.68%, digamos grosso modo 96%:

vpa(log(0.04321)/-omega_n/xi,5)

% O, puestos a decir, el transitorio dura (igualando la exponencial a 0.0183), criterio 98.2%

vpa(log(0.0183)/-omega_n/xi,3)

%% Oscilaciones: frecuencia propia, tiempo de pico, valor máximo
% A la frecuencia de las oscilaciones se le denomina frecuencia PROPIA (y es, claro, la parte IMAGINARIA de los polos):

w_propia=omega_n*sqrt(1-xi^2)

% La "velocidad" (los sistemas de 2o orden mecánicos su estado usualmente es posición y velocidad):

velo=simplify(diff(Y,t))

% Nada muy interesante que sacar, excepto que se iguala a cero (máximos y mínimos de la posición) cuando el seno es cero, o sea en

solve(velo==0,t,"ReturnConditions",true)
t_1erpico=pi/w_propia

% El valor de la salida en ese primer pico es:

Ymaxima=simplify(Y(t_1erpico))

% Al rebasamiento (porcentual) sobre el valor final se le denomina sobreoscilación (usualmente llamado \delta en los libros de texto):

SobreOsc=(Ymaxima-A*K)/(A*K)

%% Ejemplo numérico

s=tf('s');
K=3; xi=0.4; omega_n=3;
G=K*omega_n^2/(s^2+2*xi*omega_n*s+omega_n^2)

% Vamos a ver la respuesta, y luego comprobaremos fórmulas

step(G,7), grid on 

% Comprobemos todo lo dicho teóricamente:

pole(G)
sgm=xi*omega_n %parte real, bueno, cambiada de signo.
t_est_95=3/sgm
t_est_96=pi/sgm
t_est_98=4/sgm
w_propia=sqrt(1-xi^2)*omega_n %parte imaginaria
periodooscil=2*pi/w_propia
valorfin=K*1
sobreoscilacion=exp(-xi*pi/sqrt(1-xi^2))
valormax=(1+sobreoscilacion)*valorfin
t_primerpico=pi/w_propia
step(G,7), grid on %A=1   K=3, tau=5   
hold on
yline(valorfin,'-.r',Label="valor final")
yline(valormax,'-.g',Label="Valor Máximo")
yline(0.95*valorfin)
xline(t_primerpico,'-.k',Label="tiempo primer pico",LabelVerticalAlignment="bottom")
xline(t_est_95,'-.',Label="95%",Color=[.5 .5 .5],LabelVerticalAlignment="bottom")
xline(t_est_96,'-.',Label="96%",Color=[.5 .5 .5],LabelVerticalAlignment="middle")
xline(t_est_98,'-.',Label="98%",Color=[.5 .5 .5],LabelVerticalAlignment="bottom")
%Envolventes exponenciales (no te preocupes mucho de ello)
envex=1+xi^2/(1-xi^2);
fplot(3*(1-envex*exp(-sgm*t)),[0 7],Color=[.8 .9 .7],LineStyle=':',LineWidth=2)
fplot(3*(1+envex*exp(-sgm*t)),[0 7],Color=[.8 .9 .7],LineStyle=':',LineWidth=2)
hold off

%
%% Otros parámetros de posible interés (tpo. retardo / subida)
% Tiempo de "retardo": si lo hay, tiempo que tarda a "empezar" a subir... con ruido, se suele decir "a llegar al 10% del valor final".
% Tiempo de "subida": tiempo para llegar del 10 al 90% del valor final por primera vez.
% Podrían tener fórmula descrita en algún libro, pero no son comunes en textos introductorios y, al final, si hay tanto que memorizar... pues lo simulo y punto, se ve en la gráfica:

step(G,3.6), grid on
yline(0.1*valorfin,'-.')
yline(0.9*valorfin,'-.')
xline(0.15,'r')
xline(0.65,'r')

%% Conclusión del ejemplo: descripción en términos "útiles" de la respuesta
% El sistema objeto de estudio:
% Alcanza un nuevo equilibrio (estable) que incrementa salida en 3 unidades por cada unidad de incremento de entrada.
% tarda 0.15 segundos en "notarse claramente que sube" (tiempo de retardo)
% tarda 0.5 segundos más en llegar al 90% del salto final (tiempo de subida, o sea, alcanzará en 0.65 el 90%).
% tarda 1.14 segundos en alcanzar el máximo, que rebasa un 25% el valor final (tiempo de pico, sobreoscilación porcentual).
% Oscila a 2.75 rad/s durante el transitorio, período 2.3 segundos.
% El transitorio dura entre 2.6 y 3.3 segundos, dependiendo del criterio para considerarlo "terminado"... o sea, que vamos a ver un período y medio de oscilaciones y luego desaparecen, aproximadamente.
%