Transformada Laplace – Python

Transformadas de Laplace

La librería Sympy ofrece algunas opciones para revisar el tema de transformadas de laplace. Por ejemplo se obtienen los siguientes resultados a partir de funciones en el tiempo.

 f(t) = e**(2t) 

 2*t
e   
   1       s      
(-----, 2, - != 1)
 s - 2     2      

 f(t) = Dirac delta: 

DiracDelta(t - 2)
  -2*s            
(e    , -oo, True)

 f(t) = Piecewise: 

/0  for t < 0 
< \1 for t >= 0
 1          
(-, 0, True)
 s          

 f(t) = Heaviside: 

Heaviside(t)
 1          
(-, 0, True)
 s          
>>>

Las respuestas tienen tres partes (f,a,cond) : la función f(s), una constante a que define el plano de convergencia Re(s)>a y una condición auxiliar para convergencia.

usando las intrucciones, tiene una guia para continuar con las otras expresiones.

Note que para el impulso se dan las opciones de Piecewise() y Heaviside(). En las respuestas con escalones unitarios, Sympy usa Heaviside().

# Transformadas de Laplace con Sympy

import sympy as sym
t = sym.Symbol('t')
s = sym.Symbol('s')

# en dominio tiempo
e2 = sym.exp(2*t)
# escalon unitario
u = sym.Piecewise((0,t<0),(1,t>=0))
uH = sym.Heaviside(t)
# impulso
delta2 = sym.DiracDelta(t-2)

# Transformando a Laplace
e2s = sym.laplace_transform(e2,t,s)
delta2s = sym.laplace_transform(delta2,t,s)
us = sym.laplace_transform(u,t,s)
uHs = sym.laplace_transform(uH,t,s)

# Salida
print('\n f(t) = e**(2t) \n')
sym.pprint(e2)
sym.pprint(e2s)
print('\n f(t) = Dirac delta: \n')
sym.pprint(delta2)
sym.pprint(delta2s)
print('\n f(t) = Piecewise: \n')
sym.pprint(u)
sym.pprint(us)
print('\n f(t) = Heaviside: \n')
sym.pprint(uH)
sym.pprint(uHs)

Transformadas Inversas de Laplace

Se aplica el proceso contrario a las funciones anteriores y el ejemplo para la respuesta al impulso de la sección del sistema Modelo.

 H(s) = 

     s      
------------
 2          
s  + 3*s + 2

 / t    \  -2*t             
-\e  - 2/*e    *Heaviside(t)

 u(s) = 

 1          
(-, 0, True)
 s          

Heaviside(t)

 uH(s) = 

 1          
(-, 0, True)
 s          

Heaviside(t)
>>> 

En las instrucciones de  Python, donde se usan las respuestas anteriores de la transformada, se toma solo la primera parte (ejemplo:us[0]).

# Transformada inversa de Laplace
Hs = s/(s**2+3*s+2)

Ht = sym.inverse_laplace_transform(Hs,s,t)
ut = sym.inverse_laplace_transform(us[0],s,t)
uHt = sym.inverse_laplace_transform(uHs[0],s,t)

print('\n ******\n')
print('\n H(s) = \n')
sym.pprint(Hs)
sym.pprint(Ht)
print('\n u(s) = \n')
sym.pprint(us)
sym.pprint(ut)
print('\n uH(s) = \n')
sym.pprint(uHs)
sym.pprint(uHt)

Transformada Laplace

Para una señal X(t), Su transformada de Laplace esta definida como:

X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} \delta t

la señal x(t) es la inversa de la transformada X(s), que se obtiene como:

x(t) = \frac{1}{2πj} \int_{c-j\infty}^{c+j\infty} x(t) e^{st} \delta t

donde c es una constante seleccionada para asegurar la convergencia de la integral.

Los pares de ecuaciones conocidos como «Pares de la transformada de Laplace» se escriben simbólicamente como:

X(s) = L[x(t)] x(t) = L^{-1}[X(s)]

y propiedades como

L[a_1 x_1(t) + a_2 x_2 (t)] = a_1 X_1(s) + a_2 X_2 (s)

que muestra la linealidad de la transformada

Región de convergencia (ROC)

También conocida como región de existencia, es el grupo de valores de s, en la region del plano complejo, para los cuales la integral de la transformada converge.

Ejemplo

referencia: Lathi ejemplo 4.1. pdf231

Para una señal x(t) = e-atu(t), encuentre la transformada X(s) y su región de convergencia (ROC)

X(s) = \int_{-\infty}^{\infty} e^{-at} \mu (t) e^{-st} \delta t

pero tomando en cuenta que u(t)=0 para t<0 y u(t)=1 para t≥0:

X(s) = \int_{0}^{\infty} e^{-at} \mu(t) e^{-st} \delta t X(s) = \int_{0}^{\infty} e^{-at} e^{-st} \delta t = \int_{0}^{\infty} e^{-(s+a)t} \delta t = \frac{1}{s+a} \Big|_{0}^{\infty}

Si s es un número complejo y t →∞, el término exponencial no necesariamente desaparece.

El número complejo se puede expresar como:

z = α + jβ
e-zt = e-(α + jβ)t = e-αte-jβt.

El valor |e-jβt| = 1 sin importar el valor de βt.
Entonces si t→∞, e-zt →0 solo si α>0
y e-zt →∞ en el caso que α<0.

La región de convergencia (ROC) de X(s) es que la parte real s>-a.

Tabla de transformadas de Laplace (wikipedia)

Tabla de transformadas y propiedades de Laplace (Universidad de Córdova)

Propiedad de diferenciación

\frac{\delta x}{\delta t} \Leftrightarrow sX(s) - x(0^{-}) \frac{d^2x}{dt^2} \Leftrightarrow s^2X(s) - sx(0^{-}) - x'(0^{-}) \frac{\delta^3x}{\delta t^3} \Leftrightarrow s^3 X(s) - s^2 x(0^{-}) - sx'(0^{-}) - x''(0^{-})

Transformada Laplace Ejercicio01

Ejemplo 1:

Referencia: Oppenheim ejemplo 9.37 p718/pdf746, semejante al ejercicio sistema-modelo LTI Lathi 1.8-1 pdf/p.80. Oppenheim problema 2.61c pdf/p.191

Suponga un sistema LTI causal que se describe mediante la ecuación diferencial:

\frac{d^2y(t)}{dt} + 3 \frac{dy(t)}{dt} + 2y(t) = x(t)

junto con la condición de reposo inicial:

H(s) = \frac{1}{s^2 + 3s +2}

Suponga que x(t) = α μ(t) , entonces la entrada del sistema es x(s) = α

Siempre que las señales de entrada son cero para t<0, la convolución es la multiplicación de los términos. (no aplica para valores t<0)

y(s) = H(s) X(s) = \Big[\frac{1}{s^2 + 3s +2} \Big]\Big[\frac{\alpha}{s}\Big] y(s) = \frac{\alpha}{2}\frac{1}{s} - \alpha\frac{1}{s+1} +\frac{\alpha}{2}\frac{1}{s+2}

por lo que aplicando la anti transformada:

y(t) = \alpha \Big [ \frac{1}{2} - e^{-t} + \frac{1}{2} e^{- 2t} \Big] u(t)

Usando Python y Sympy

Para crear la expresión polinómica simple, no usamos el factor α, que se añade luego de obtener el resultado, se obtiene:

    a         a      a 
--------- - ----- + ---
2*(s + 2)   s + 1   2*s

resultado al que hay que multiplicar por la constante α todos sus términos.

las instrucciones son:

# Oppenheim ejemplo 9.37 p718/pdf746
import sympy as sym

s = sym.Symbol('s')
a = sym.Symbol('a')
Hs = (1/(s**2+3*s+2))*(a/s)
separa = Hs.apart(s)

sym.pprint(separa)

Ejemplo 2

Referencia: Oppenheim ejemplo 9.38 p719/747

Considere el sistema caracterizado por la ecuación diferencial con condiciones iniciales.

\frac{d^2y(t)}{dt^2} + 3 \frac{dy(t)}{dt} + 2y(t) = x(t) y(0^{-}) = \beta \text{, } y'(0^{-}) = \gamma

sea x(t) = αμ(t)

Aplicando la transformada de laplace en ambos lados de la ecuación con condiciones iniciales y aplicando la propiedad de diferenciación:

\frac{d^2y(t)}{dt^2} = s^2Y(s) - y(0^{-})s - y'(0^{-}) = s^2Y(s) - \beta s - \gamma \frac{dy(t)}{dt} = sY(s) - \beta

resultados que se insertan en la ecuacion original:

[s^{2}Y(s) - \beta s - \gamma] + 3[sY(s) - \beta] + 2Y(s) = \frac{\alpha}{s}

para simplificar como:

[s^{2}+3s +2] Y(s) - \beta s - (\gamma + 3 \beta) = \frac{\alpha}{s} [s^{2}+3s +2] Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) (s+1)(s+2)Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) Y(s) = \frac{\alpha}{s(s+1)(s+2)}+ \frac{\beta s}{(s+1)(s+2)} + \frac{(\gamma + 3 \beta)}{(s+1)(s+2)}

se encuentra que:

Y(s) = \frac{\beta (s+3)}{(s+1)(s+2)} + \frac{\gamma}{(s+2)(s+2)} + \frac{\alpha}{s(s+2)(s+2)}

donde Y(s) es la transformada unilateral de Laplace de y(t)

el último término, representa la respuesta del sistema LTI causal y la condición de reposo inicial, conocido también como la respuesta en estado cero.

Una interpretación análoga se aplica a los primeros dos términos del miembro derecho de la ecuación, que representa la respuesta del sistema cuando la entrada es cero (α=0) , conocida también como respuesta a entrada cero.

Observe que la respuesta a entrada cero es una función lineal de los valores de las condiciones iniciales. Demostrando que la respuesta es la superposición del estado cero y respuesta a entrada cero.

si α=2, β = 3 y γ=-5, al realizar la expansión en fracciones parciales encontramos que:

Y(s) = \frac{1}{s} - \frac{1}{s+1} + \frac{3}{s+2}

con la aplicación de las tablas de transformadas se obtiene:

y(t) = \Big[ 1 - e^{-t} + 3 e^{-2t} \Big] \mu (t) \text{, para } t\gt 0

usando Python y Sympy

se puede trabajar la ecuación para obtener:

ecuacion de entrada
    2                                   a
Ys*s  + 3*Ys*s + 2*Ys - b*s - 3*b - g = -
                                        s
ecuacion Ys:  
        2               
 a + b*s  + 3*b*s + g*s 
[----------------------]
      / 2          \    
    s*\s  + 3*s + 2/    
fracciones parciales de Ys: 
 a    a - 2*b - 2*g   a - 2*b - g
--- + ------------- - -----------
2*s     2*(s + 2)        s + 1   
sustituye a=2, b= 3 y g=-5
  3       1     1
----- - ----- + -
s + 2   s + 1   s
>>>

con las siguientes instrucciones, presentadas por cada sección de la parte de desarrollo analítico, para ir comparando los resultados.

# Oppenheim ejemplo 9.38 p719/747
import sympy as sym

s,a,b,g  = sym.symbols('s a b g')
Ys = sym.Symbol('Ys')

d2y = (s**2)*Ys - b*s - g
dy  = s*Ys - b

ecuacion = sym.Eq(d2y + 3*dy + 2*Ys, a/s)
ecuacionYs = sym.solve(ecuacion,Ys)
parciales = ecuacionYs[0].apart(s)

solucion = parciales.subs({a:2,b:3,g:-5})

print('ecuacion de entrada')
sym.pprint(ecuacion)
print('ecuacion Ys:  ')
sym.pprint(ecuacionYs)
print('fracciones parciales de Ys: ')
sym.pprint(parciales)
print('sustituye a=2, b= 3 y g=-5')
sym.pprint(solucion)

La definición de la ecuacion se realiza usando las variables como símbolos, incluyendo Y(s) como Ys. En éste ejercicio no existe X(s), por lo que no se la incluyó.

Luego se busca la expresión solución para Ys, por simplicidad se la cambia a fracciones parciales, asi es más sencillo usar la tabla de transformadas de Laplace.

Para una solución con los valores de a,b y g, se usa un diccionario {} con las parejas de variables y sus valores. Compare los resultados con lo expuesto en la parte teórica.

Fracciones Parciales – Separar con Python

Referencia: Oppenheim 9.4 p671/pdf699. Lathi B5 pdf21, Schaum/HSu 3.8.C p120

Se pueden separar en fraciones parciales del siguiente ejemplo:

H(s) = \frac{2s^2 +5}{s^2 +3s+2}

se requiere obtener los residuos, polos y términos directos al separar en fracciones parciales entre dos polinomios P(s)/Q(s).

\frac{P(s)}{Q(s)} = \frac{ceros(1)}{s - polos(1)}+\frac{ceros(2)}{s - polos(2)}+\text{ ... } \text{ ... }+\frac{ceros(n)}{s - polos(n)}+ganancia(s)

Se puede realizar de dos formas usando Python, en ambas se requiere definir los vectores P y Q.

Sympy – fórmulas simbólicas

La primera usando Sympy para crear las fórmulas simbólicas y luego formar una  fracción, separar las partes Sympy.apart()

Para los datos proporcionados en el ejemplo se tiene como resultado:

P = [2, 0, 5]
Q = [1, 3, 2]

Fracciones parciales:
2 - 13/(s + 2) + 7/(s + 1)

con las instrucciones sympy:

# Expansión en fracciones parciales
# P es numerador, Q es denominador

P = [2, 0, 5]
Q = [1, 3, 2]

import sympy as sym
s = sym.Symbol('s')
n = len(P)
m = len(Q)
Ps = 0
for i in range(0,n,1):
    Ps = Ps + P[i]*s**(n-i)
Qs = 0
for i in range(0,m,1):
    Qs = Qs + Q[i]*s**(m-i)
fraccion = Ps/Qs
parciales = sym.apart(fraccion)

# Salida - polinomio
print('Fracciones parciales:')
print(parciales)

Las fórmulas simbólicas permiten agrupar y realizar operaciones entre varios sistemas además de expandir, separar, y evaluar las expresiones. Util principalmente para la sección de transformadas de Laplace.

Un resumen de sympy se muestra en:

Fórmulas simbólicas – Sympy


Scipy – Procesar señales

La segunda forma es usando libreria scipy.signal.residue con los valores de P y Q:

ceros(n)   :  [-13.   7.]
polos(n)   :  [-2. -1.]
ganancia(s):  [ 2.]
>>>

Con los resultados el equivalente en Laplace es:

transformado al dominio del tiempo,
H(t) = (−13e−2t + 7e−t)u(t) + 2δ(t)

las instrucciones adicionales a la sección anterior son:

# Usando scipy.signal.residue

import scipy.signal as senal
ceros,polos,ganancia = senal.residue(P,Q)

# SALIDA - coeficientes
print()
print('ceros(n)   : ', ceros)
print('polos(n)   : ', polos)
print('ganancia(s): ', ganancia)

Diagramas de Bloques con s

Representaciones en diagramas de bloques para los sistemas LTI causales descritos por ecuaciones diferenciales y funciones racionales del sistema. (Oppenheim 9.8.2)

Ejemplo 1.

referencia: Oppenheim ejemplo 9.28 p708/pdf736

Considere un sistema LTI causal

\frac{\delta y(t)}{\delta t} + 3y(t) = x(t) Dy(t) + 3y(t) = x(t) Dy(t) = x(t) -3y(t) y(t) = \frac{1}{D}[x(t) - 3y(t)]

usando s en lugar de D

y(s) = \frac{1}{s} [x(s) - 3 y(s)]

usando la transformada de Laplace, cambiando s por d/dt desde la primera ecuación, la función del sistema es:

( s+3) y(s) = x(s)

por lo que y(s)/x(s) es:

H(s) = \frac{1}{s+3}

al cambiar el signo en el coeficiente de la multiplicación de la trayectoria de retroalimentación, obtenemos una forma de diagrama de bloques similar a la anterior, pero con la que se verifica que:

H(s) = \frac{\frac{1}{s}}{1+\frac{3}{s}} = \frac{1}{s+3}

que confirma que:

\frac{Y(s)}{X(s)} = H(s) = \frac{H_1(s)}{1+H_1(s)H_2(s)}

siendo H1 el bloque superior, y H2 el bloque inferior en el diagrama.


Ejemplo 2

referencia: Oppenheim ejemplo 9.30 p711/pdf739

Considere el sistema causal de segundo orden con la funcion del sistema:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{1}{s^2 + 3s +2}

la entrada x(t) y la salida y(t) para este sistema satisfacen la ecuación diferencial:

\frac{\delta^2(t)}{\delta t^2} + 3 \frac{\delta y(t)}{dt} + 2 y(t) = x(t) (s^2 + 3s + 2)y(t) = x(t)

pero también se puede armar como:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{k_1}{s+1} + \frac{k_2}{s+2}

siendo el numerador :

k_1s + 2k_1 + k_2s + k_2 = (k_1 + k_2)s + (2k_1 +k_2)

que se observa que debe dar:
(k_1 + k_2)s + (2k_1 +k_2) = 0s + 1

entonces

k_1 + k_2 = 0 2k_1 + k_2 = 1

con lo que

k1 = 1 y
K2 = -1
H(s) = \frac{1}{s+1} -\frac{1}{s+2}

Si

H(s) = \frac{1}{(s+1)(s+2)}

se escribe como

H(s) = \Big[\frac{1}{s+1}\Big]\Big[\frac{1}{s+2}\Big]

que se representa como la serie o cascada de dos sistemas de primer orden.

que por cierto, también es la convolución de dos sistemas en serie.

También es posible realizar las fracciones parciales con Sympy y se obtiene:

    1       1  
- ----- + -----
  s + 2   s + 1
>>>

usando las instrucciones:

# Oppenheim 9.30 p711/pdf739
import sympy as sym

s = sym.Symbol('s')
Hs = 1/((s+1)*(s+2))
separa = Hs.apart()

sym.pprint(separa)

Ejemplo 3.

referencia: Oppenheim ejemplo 9.31.p712/pdf739

Considere la funciín del sistema

H(s) = \frac{2s^2 + 4s - 6}{s^2 + 3s + 2}

se puede reescribir como:

H(s) = \Big[\frac{1}{s^2 + 3s + 2} \Big] \Big[ 2s^2 + 4S - 6 \Big]

y si las derivadas se toman del sistema anterior, se tiene una forma mas simplificada:

también, con un poco de trabajo, reordenando las ecuaciones, se podría convertir en:

H(s) = \Big[ 2\frac{s-1}{s+2} \Big] \Big[\frac{s+3}{s+1} \Big]

o también en :

H(s) = 2 + \frac{6}{s+2} - \frac{8}{s+1}

que requiere un sistema en forma paralelo.
QUEDA COMO TAREA HACER EL DIAGRAMA.

Se presenta el ejercicio usando Sympy:

      6       8  
2 + ----- - -----
    s + 2   s + 1
>>>

que se obtiene usando las instrucciones:

# Oppenheim 9.31.p712/pdf739
import sympy as sym

s = sym.Symbol('s')
Hs = (2*s**2+4*s-6)/(s**2+3*s+2)
separa = Hs.apart()

sym.pprint(separa)