%% Respuesta temporal: fracciones simples y transformada inversa de Laplace con retardos,
%% a) sistema 1er orden con retardo
% Presentaciones vídeo: http://personales.upv.es/asala/YT/V/dly1er1.html , http://personales.upv.es/asala/YT/V/dly1er2.html
%% b) sistema 1er orden inestable ante pulso de entrada
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/dly1erpul.html
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Este código ejecutó sin errores en Matlab R2022b (Linux)
% Objetivo: revisar diferentes ejemplos de cálculo de respuesta temporal por transformada de Laplace inversa usando la descomposición en fracciones simples (partial fraction expansion) de la Symbolic Math Toolbox de Matlab. En este caso, con presencia de retardos en función de transferencia, entrada o ambas.
%
% Tabla de Contenidos
% Respuesta temporal: fracciones simples y transformada inversa de Laplace con retardos,
% a) sistema 1er orden con retardo
% b) sistema 1er orden inestable ante pulso de entrada
%     Ejemplo 1: Sistema de 1er orden con retardo ante escalón, c.i. = 0
%         Método 1: solucionar sub-problema sin retardo y luego retrasar la solución
%         Método 2: por Laplace manejando exponenciales en "s":
%             Descomposición en fracciones simples "manual"
%             Descomposición en fracciones simples "automática"
%             Transformada Inversa de Laplace
%     Ejemplo 2: Sistema inestable 1er orden ante pulso rectangular, c.i. no nulas
%         Método 1 (intuitivo y posiblemente lo que yo recomiende): solucionar "fragmentos" sin retardo y luego ya retrasar y superponer.
%         Método 2: Laplace con operador retardo (exp(-ds)).
%             Descomposición en fracciones simples "manual"
%             Descomposición en fracciones simples "automática"

syms y(t) %función del tiempo simbólica
syms s %variable de Laplace
syms A B real %incógnitas varias en descomposiciones "a mano"

%
%% Ejemplo 1: Sistema de 1er orden con retardo ante escalón, c.i. = 0
% Consideremos:     \frac{dy}{dt}+3y(t)=u(t-2) ante entrada escalón de 6 unidades, suponiendo y(0)=0.
% Usemos "dsolve" para ver qué tiene que dar al final. El escalón en symbolic toolbox es "heaviside"; es necesario usarlo cuando no está aplicado en t=0.

sol=dsolve(diff(y)+3*y==6*heaviside(t-2), y(0)==0); %esto debe dar al final.
simplify(sol,100)
fplot(sol,[0 5],LineWidth=3), grid on

%% Método 1: solucionar sub-problema sin retardo y luego retrasar la solución
% Básicamente es la solución sin retardo, recorrida 2 unidades a la derecha... esto se entiende un poco mejor, quizás:

sol_sinretardo=dsolve(diff(y)+3*y==6, y(0)==0) 

% Será, al retrasar
% y(t)=\begin{cases} 0 &t<2\\ 2-2e^{-3(t-2)} &t\geq 2\end{cases}
% Podríamos resolver y'+3y=6 por Laplace (o lo que sea) y luego "retrasarlo". Posiblemente, ¡hasta lo recomiendo!.
% No lo detallo, me creo "dsolve" porque problemas similares se han resuelto en otros materiales. Con Matlab retrasaríamos como:

yretrasada=subs(heaviside(t)*sol_sinretardo,t,t-2)

% Gráficamente:

fplot(sol_sinretardo,[0 5],LineWidth=2,LineStyle="-.",Color=[.5 .6 .7]), grid on, hold on
fplot(yretrasada,[0 5],LineWidth=3), hold off %coincide con dsolve, obviamente
legend("solución sin retardo","solución retardada correcta",Location="best")

%
%
%% Método 2: por Laplace manejando exponenciales en "s": e^{-ds}
% Por ser "puristas" siguiendo las metodologías en muchos libros, vamos a hacerlo por Laplace con "la transformada de Laplace del operador retardo", que es e^{-ds} siendo d\geq 0 el tiempo de retardo. En efecto, una función f(t) "retrasada" se escribiría
% f_{retra}(t)=\begin{cases}0 &t<d \\ f(t-d) &t\geq d\end{cases}
% de modo que, usando el cambio [t-d=\tau,\ t=\tau+d]:
% \mathcal L(f_{retra})=\int_d^\infty f(t-d)e^{-st}\,dt=\int_0^{\infty} f(\tau)e^{-s\cdot (\tau+d)}\,d\tau=e^{-d\cdot s}\int_{0}^\infty f(\tau)e^{-s\tau}\,d\tau=e^{-d\cdot s}\cdot \mathcal L(f)
%
% Esto de llenar de "exponenciales en s" las expresiones (ojo, no en TIEMPO!)  ya no tiene tanta interpretación "intuitiva" como el primer método, aunque es "común" entre el profesorado que explica este tema (yo mismo, je)... Bueno, al final, si entiendes que las dos cosas son esencialmente lo mismo, no es tan difícil, pero la primera vez asusta un poco y puedes confundirte entre las exponenciales "en s" y "en t".
% La transformada de Laplace de \frac{dy}{dt}+3y(t)=u(t-2) es:
% sY(s)-y(0)+3Y(s)=e^{-2s}\cdot U(s),
% esto es:
% Y(s)=\frac{e^{-2s}}{s+3}U(s)+\frac{y(0)}{s+3}

G=exp(-2*s)/(s+3); %Funcion de transferencia
U=6/s; %Transformada de Laplace de la entrada u(t)=6, constante (escalón).
%suponemos y(0)=0, por lo que no consideramos término de cond. iniciales
Y=G*U %Transformada de Laplace de la salida ante c.i.=0

%% Descomposición en fracciones simples "manual"
% Para hacerlo "manualmente"; deberíamos igualar Y(s)=\frac{A}{s}+\frac{B}{s+3}

Desco=A/s+B/(s+3)

% Pasando a común denominador:
% \frac{6e^{-2s}}{s(s+3)}=Y(s)=\frac{A}{s}+\frac{B}{s+3}=\frac{A(s+3)+Bs}{s(s+3)}

Desco2=collect(simplifyFraction(Desco),s)

% Deberíamos igualar numeradores:

sol=solve([6*exp(-2*s) 0]==[3*A  A+B],{A,B})

%
%% Descomposición en fracciones simples "automática"

partfrac(Y) %igual que arriba, claro.

%% Transformada Inversa de Laplace
% De hecho, la Symbolic toolbox hace la transformada inversa de todo (como él quiera, usando o no partfrac internamente) sin necesidad de lo de arriba:

Y_t=simplify(ilaplace(Y),30)
fplot(Y_t,[0 5],LineWidth=2), grid on,title("Salida")

%%
%% Ejemplo 2: Sistema inestable 1er orden ante pulso rectangular, c.i. no nulas
% El modelo EDO del sistema es:          \frac{dy}{dt}-\frac{1}{3}y(t)=u(t)
% Su transformada es:
% sY(s)-y(0)-Y(s)/3=U(s),
% esto es:
% Y(s)=\frac{1}{s-1/3}U(s)+\frac{y(0)}{s-1/3}

G=1/(s-1/3); %Función de transferencia

% La entrada va a ser un pulso rectangular de amplitud 0.25 y duración 2 segundos:

U_time=0.25*(heaviside(t)-heaviside(t-2))
fplot(U_time,[-2 5],LineWidth=3), grid on, xline(0,'-.r',LineWidth=2)
title("Entrada tipo pulso (vale 0 \forall t\leq 0)")
y0=-0.35; %condición inicial
sol=simplify( dsolve(diff(y)-y/3==U_time, y(0)==y0) , 50) %esto debe dar al final.

% No lo entiende "ni Dios", pero vamos a dibujarlo:

fplot(sol,[0 7],LineWidth=3,Color=[.8 .3 .1]), grid on, title("Salida")

%% Método 1 (intuitivo y posiblemente lo que yo recomiende): solucionar "fragmentos" sin retardo y luego ya retrasar y superponer.
% La entrada es la superposición de 2 escalones:

U1=0.25/s;
U2=-0.25/s*exp(-2*s); %el mismo retrasado 2 segundos, deshaciendo la subida.
fplot(ilaplace([U1+U2]),[0 5], LineWidth=4),
hold on
fplot(ilaplace([U1;U2]),[0 5], LineWidth=3, LineStyle='-.'),
hold off
grid on, legend("u_1+u_2","u_1","u_2")

% La solución ante u_1(t)=0.25, U_1(s)=0.25/s será:

Y1_s=G*U1

% Con las adecuadas fracciones simples y tablas, resultará en:

Y1_t=ilaplace(Y1_s)

% La salida ante u_2(t)=\begin{cases}0 &t< 2 \\ -0.25 &t\geq 2\end{cases}, U_2(s)=-\frac{0.25}{s}\cdot e^{-2s} será la misma que ante u_1 pero retrasada, o sea
% y_2(t)=\begin{cases}0 &t< 2 \\ -y_1(t-2) &t\geq 2\end{cases}

Y2_t=subs(-heaviside(t)*Y1_t,t,t-2)

% y, por último la respuesta libre (salida ante entrada cero y condiciones iniciales no nulas) será:

Yfree=ilaplace(y0/(s-1/3))

% con lo cual, la salida será y(t)=y_1(t)+y_2(t)+y_{free}(t), esto es:
% y(t)=-\frac{7}{20}e^{t/3}+\frac{3}{4}(e^{t/3}-1)+\begin{cases}0 &t< 2 \\ -\frac{3}{4}(e^{(t-2)/3}-1) &t\geq 2\end{cases}
% que, reorganizando queda:
% y(t)=\begin{cases}-\frac{7}{20}e^{t/3}+\frac{3}{4}(e^{t/3}-1) &t< 2 \\ -\frac{7}{20}e^{t/3}+\frac{3}{4}(e^{t/3}-1)-\frac{3}{4}(e^{(t-2)/3}-1) &t\geq 2\end{cases}
% *Usando Matlab,

Y_metodo1=simplify(Yfree+Y1_t+Y2_t)

% coincide con dsolve, como era de esperar.

fplot(Y_metodo1,[0 7]), grid on

%
%% Método 2: Laplace con operador retardo (exp(-ds)).
%
% Ahora por Laplace, poco a poco a mano:

U=laplace(U_time)
Y=simplify(G*U+y0/(s-1/3),70)

%% Descomposición en fracciones simples "manual"
% Para hacerlo "manualmente"; deberíamos igualar Y_{forz}(s)=\frac{A}{s}+\frac{B}{s-\frac{1}{3}}

Desco = A/s +  B/(s-1/3)

% Pasando a común denominador:
% Y_{forz}(s)=\frac{A}{s}+\frac{B}{s-\frac{1}{3}}=\frac{A(s-\frac{1}{3})+Bs}{s(s-\frac{1}{3})}

Desco2=collect(simplifyFraction(Desco),s)

% Como los denominadores coinciden el de Yforz y el de Desco2, deberíamos igualar numeradores:

sol=solve([15/20*(1-exp(-2*s)) -21/20]==[-A  3*(A+B)],{A,B})

%
%% Descomposición en fracciones simples "automática"
% La Symbolic Toolbox lo hace en una línea, lo de arriba "sobra" si tenemos Matlab a mano...

partfrac(Y)
Y_t=ilaplace(Y)
fplot(Y_t,[0 7]), grid on, yline(0)