%% Análisis de la respuesta temporal: régimen estacionario ante rampa unitaria de sistema, dada su función de transferencia
% Presentaciones en vídeo:
%  --- http://personales.upv.es/asala/YT/V/regestramp1.html
%  --- http://personales.upv.es/asala/YT/V/regestramp2.html
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Este código funcionó sin errores en Matlab R2022b (Linux)
% Objetivo: calcular el régimen estacionario P+Qt de un sistema con FdT G(s) ante la entrada u(t)=t.
% Alternativamente, la interpretación como "traslación temporal" Q\cdot (t+\delta), con \delta=P/Q.
% Obviamente, ante una entrada u(t)=At todo se multiplica por A, por linealidad, por eso sólo estudiamos rampa "unitaria" A=1.
%
% Tabla de Contenidos
% Desarrollo teórico
%     A partir de la descomposición en fracciones simples
%     Deducción alternativa (serie de Taylor)
%     Interpretación como retardo/adelanto temporal
% Aplicación a sistemas de primer y segundo orden
%     Sistemas de primer orden
%     Sistemas de segundo orden
% Ejemplos Numéricos
%     Ejemplo 1
%     Ejemplo 2
%     Ejemplo 3
%
%% Desarrollo teórico
%% A partir de la descomposición en fracciones simples
%
% Consideremos G(s)=\frac{N(s)}{D(s)}, sin polos en el origen.
% La salida ante rampa unitaria en Laplace es:
% y(s)=\frac{N(s)}{D(s)}\cdot \frac{1}{s^2} = \frac{P}{s}+\frac{Q}{s^2}+\frac{R(s)}{D(s)}=\frac{P\cdot sD(s)+Q\cdot D(s)+s^2R(s)}{D(s)s^2}
% se verificará:
% N(s)=P\cdot sD(s)+Q\cdot D(s)+s^2R(s)
% Y por tanto, haciendo s=0 tendremos:
% N(0)=QD(0), o sea Q=G(0).
% Por otra parte, haciendo derivadas respecto a s, la igualdad debería seguir siendo cierta:
% N'(s)=P\cdot D(s)+P\cdot sD'(s)+Q\cdot D'(s)+s^2R'(s)+2sR(s)
% de modo que haciendo, de nuevo s=0, daría lugar a:
% N'(0)=P\cdot D(0)+Q\cdot D'(0)
% Sustituyendo el valor de Q=G(0)=N(0)/D(0) calculado antes, tenemos:
% N'(0)=P\cdot D(0)+Q\cdot D'(0) = P\cdot D(0)+\left(\, N(0)/D(0)\,\right) \cdot D'(0)
% O sea:
% N'(0)\cdot  D(0) = P\cdot D(0)^2+ N(0) \cdot D'(0)
% resultando en:
% P=\frac{N'(0)\cdot D(0)-N(0)\cdot D'(0)}{D(0)^2}=G'(0)
% De modo que aparte del transitorio con las mismas exponenciales que la respuesta libre R(s)/D(s), habrá un régimen estacionario:
% y_{estac}(t)=P+Qt = \frac{dG}{ds}(0)+G(0)\cdot t
%
%
%
%% Deducción alternativa (serie de Taylor)
% Si G(s) no tiene polos en el origen, haremos Taylor:
% G(s)=G(0)+G'(0)s+X(s)s^2
% Y(s)=G(s)/s^2=G(0)/s^2+G'(0)/s+X(s)
% y(t)=G(0)\cdot t + G'(0) \ + \ \text{[términos exponenciales]}
% Con lo que, identificando X(s)=\frac{R(s)}{D(s)} tenemos una segunda deducción del régimen permanente ante rampa.
%
%
%
%% Interpretación como retardo/adelanto temporal
% Escribamos G(s)=K\cdot G_1(s), siendo G_1(s) tal que tiene ganancia estática unidad: G_1(0)=1, o sea que definimos K:=G(0).
% La respuesta estacionaria ante rampa de G_1(s) será G_1(0)t+G_1'(0)= t+G_1'(0), o sea que es la rampa unitaria de entrada adelantada G'_1(0) o retrasada si G'_1(0)<0.
% Por tanto, la respuesta rampa estacionaria del sistema original puede ser interpretada como una "ganancia proporcional K" y un desplazamiento temporal \delta =G_1'(0):
% y(t)=K\cdot (t+\delta)\ +  \ \text{[términos exponenciales]}
%
%% Aplicación a sistemas de primer y segundo orden
%% Sistemas de primer orden

syms K tau s
G1(s)=1/(tau*s+1)
G(s)=K*G1(s)
dG1=diff(G1,s)
dG1(0)

% La respuesta estacionaria ante rampa será K\cdot (t-\tau) la proporcionalidad K retrasada \tau segundos.
%% Sistemas de segundo orden

syms xi omega_n
G1(s)=omega_n^2/(s^2+2*xi*omega_n*s+omega_n^2)
G(s)=K*G1(s)
dG1=simplify(diff(G1,s))
dG1(0)

% La respuesta rampa estacionaria será la proporcionalidad K , retrasada \frac{2\xi}{\omega_n}:
%  y(t)=K\cdot (t-\frac{2\xi}{\omega_n})\ + \ \text{[términos exponenciales]}
%% Ejemplos Numéricos
%% Ejemplo 1

G(s)=5/(2*s+1);
K=G(0)
G1(s)=G(s)/K
dG1ds=diff(G1,s);
delta=dG1ds(0)
fplot([ilaplace(G/s^2);K*(t+delta);K*t],[0 10]), grid on

%% Ejemplo 2

G(s)=4/(s^2+1.5*s+9)
K=G(0)
omega_n=3;
xi=1.5/2/omega_n
G1=G/K;
dG1=simplify(diff(G1,s))
delta=dG1(0)
-2*xi/omega_n %coincide

% La respuesta está retrasada ( \delta<0 ) respecto a la pura ganancia (proporcionalidad instantánea). El tiempo de establecimiento (aprox.) será:

t_est=pi/(xi*omega_n)
syms t real
fplot([ilaplace(G/s^2); K*(t+delta); K*t ],[0 6]), grid on
xline(t_est,':')

%% Ejemplo 3

G(s)=(2+3*s)/(s^3+2*s^2+s+1)/(0.05*s+1)
[Num,Den]=numden(G);
vpa(solve(Den==0),4) %estable

% El régimen permanente se alcanzará aproximadamente, tras:

t_est=pi/0.1226 %parte real polos dominantes

% Pasemos a calcularlo:

K=G(0)
G1=G/K
dG1=simplify(diff(G1,s))
delta=dG1(0)

% En este caso, la respuesta estacionaria está "adelantada" respecto a la proporcionalidad instantánea (ganancia estática)... los ceros de fase mínima "adelantan"...

syms t real
fplot([ilaplace(G/s^2); K*(t+delta); K*t ],[0 35]), grid on
xline(t_est,':')

%
% *Nótese que ante "rampa" las fórmulas de tiempo de establecimiento parecen excesivamente "pesimistas" porque se ve en las gráficas que las cosas acaban antes... Claro, la rampa "sube" y por tanto, las desviaciones iniciales de la exponencial son cada vez más pequeñas... Pero no vamos a discutir otra fórmula para el tiempo de establecimiento... Las exponenciales tardan lo mismo a perder el 95% ante escalón, rampa, o libre... otra cosa es cómo de pequeñas sean ante la solución estacionaria (constante, rampa, senoidal). Ante la duda, se simula y se sacan conclusiones sobre la gráfica.