Transformada de Laplace – Tabla de Propiedades

Referencia: Oppenheim tabla 9.3 p717 pdf 748. Lathi tabla 4.2 p361. Hsu tabla 3-2 p119

Propiedad Señal Transformada de Laplace
x(t)
x1(t)
x2(t)
X(s)
X1(s)
X2(s)
Linealidad ax1(t) + bx2(t) a X1(s) + b X2(s)
Desplazamiento en t x(t-t0) e-st0X(s)
Desplazamiento
en domino de s
es0tx(t) X(s-s0)
Escalamiento
en tiempo
x(at), a>0 \frac{1}{a} X\Big(\frac{s}{a}\Big)
inversión
en tiempo
x(-t) X(-s)
Conjugación x*(t) X*(s)
Convolución
(suponiendo que
x1(t) y x2(t) son cero
para t<0)
x_1 (t) \circledast x_2 (t) X1(s)X2(s)
Diferenciación en
el dominio t
\frac{\delta}{\delta t} x(t) sX(s)-x(0)
\frac{\delta^2}{\delta t^2} x(t) s2X(s)-sx(0)-x(0)
\frac{\delta^3}{\delta t^3} x(t) s3X(s) -s2x(0)-sx(0)-x'(0)
\frac{\delta^n}{\delta t^n} x(t) s^n X(s) -\sum_{k=1}^{n}s^{(n-k)}x^{(k-1)} (0^{-})
Diferenciación en
el dominio s
-tx(t) \frac{\delta}{\delta s} X(s)
Integración en
en el dominio t
\int_{0^-}^t x(\tau) \delta \tau \frac{1}{s} X(s)
\int_{-\infty}^t x(\tau) \delta \tau \frac{1}{s} X(s)+\frac{1}{s} \int_{-\infty}^{0^{-}} x(t) \delta t
Convolución en frecuencia x1(t) x2(t) \frac{1}{2 \pi j} X_1(s) \circledast X_2(s)
valor inicial x(0+) \lim _{s \rightarrow \infty} s X(s)
n>m
valor final x(∞) \lim _{s \rightarrow 0} s X(s)
polos de sX(s) en LHP (izquierda)

Transformada de Laplace – Tabla

Transformada de Laplace – Tabla

Referencia: Lathi Tabla 4.1 p334. Oppenheim Tabla 9.2 p692. Hsu Tabla 3-1 p115

No. x(t) X(s) ROC
1a δ(t) 1 Toda s
1b δ(t-T) e-sT Toda s
2a μ(t) \frac{1}{s} Re{s}>0
2b -μ(-t) \frac{1}{s} Re{s}<0
3 tμ(t) \frac{1}{s^2} Re{s}>0
4a tnμ(t) \frac{n!}{s^{n+1}} Re{s}>0
4b \frac{t^{n-1}}{(n-1)!} \mu (t) \frac{1}{s^n} Re{s}>0
4c -\frac{t^{n-1}}{(n-1)!} \mu (-t) \frac{1}{s^n} Re{s}<0
5 eλtμ(t) \frac{1}{s-\lambda} Re{s}>0
6 teλtμ(t) \frac{1}{(s-\lambda)^2} Re{s}>0
7 tneλtμ(t) \frac{n!}{(s-\lambda)^{n+1}}
8a cos (bt) μ(t) \frac{s}{s^2+b^2} Re{s}>0
8b sin (bt) μ(t) \frac{b}{s^2+b^2} Re{s}>0
9a e-atcos (bt) μ(t) \frac{s+a}{(s+a)^2+b^2} Re{s}>-a
9b e-atsin (bt) μ(t) \frac{b}{(s+a)^2+b^2} Re{s}>-a
10 \mu_n (t) = \frac{\delta ^n}{\delta t^n} \delta (t) sn Toda s
11 \mu_{-n} (t) = \mu (t) \circledast \text{...} \circledast \mu (t)

n veces

\frac{1}{s^n} Re{s}>0
12a re-atcos (bt+θ) μ(t) \frac{(r\cos (\theta)s + (ar \cos (\theta) - br \sin (\theta))}{s^2+2as+(a^2+b^2)}
12b re-atcos (bt+θ) μ(t) \frac{0.5 re^{j \theta}}{s+a-jb} + \frac{0.5 re^{-j \theta}}{s+a+jb}
12c re-atcos (bt+θ) μ(t) \frac{As+B}{s^2+2as+c}
r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} \theta = \tan ^{-1} \Big( \frac{Aa-B}{A\sqrt{c-a^2}}\Big)

b = \sqrt{c-a^2}

12d e^{-at}\Bigg[A \cos (bt) + \frac{B-Aa}{b} \sin (bt) \Bigg] \mu (t) \frac{As+B}{s^2 + 2as+c}
b = \sqrt{c-a^2}

Transformada Laplace – Tabla de Propiedades

Transformada de Laplace – Concepto con Python

4.4.5 LTI Laplace – Circuitos electrónicos Activos «Op Amp»

Referencia: Lathi 4.4-1 p382,  Oppenheim 11.50 p896

Se amplian los conceptos de circuitos pasivos analizados con transformadas e Laplace, aplicando a circuitos activos. Se obtienen los circuitos equivalentes o modelos matemáticos y se repiten los procedimientos anteriores.

El elemento activo más conocido es el amplificador operacional  (op amp) que tienen ganancia «muy alta». El  voltaje de salida v2 = – A v1, donde A es del orden de 105 o 106. Un factor importante es que la impedancia de entrada es muy alta del orden de 1012Ω y la impedancia de salida es muy baja (50-100Ω)

La configuración de la ganancia se establece con los resistores Ra y Rb y la forma de conectar las entradas y salidas.

K = 1+\frac{R_a}{R_b} v_2 = K v_1 v_2 = (R_b + R_a) i_o = R_b i_o + R_a i_o v_2 = v_s = Ra i_o = R_a i_o \frac{v_2}{v_1} =\frac{R_b+R_a}{R_a} = 1+\frac{R_b}{R_a} = K

1. Amplificador Operacional en el dominio s

Referencia: Lathi 4.6-5 p399

Dada la alta impedancia del op amp, la corriente de retroalimentación  I(s) fluye solo por los resistores. El voltaje de entrada es cero o muy pequeño dada la ganancia muy grande del op amp. Dadas estas simplificaciones, se aproxima con mucha precisión que:

Y(s) = - I(s) Z_f(s) I(s) = \frac{X(s)}{Z(s)} Y(s) = -\frac{Z_f(s)}{Z(s)} H(s) = -\frac{Z_f(s)}{Z(s)}

1,1 Amplificador Operacional como multiplicador escalar

H(s) = -\frac{R_f}{R}

1,2 Amplificador Operacional como Integrador

Referencia: Oppenheim 11.52 p898

H(s) = \Big(-\frac{1}{RC}\Big) \frac{1}{s}

1,3 Amplificador Operacional como Sumador

Y(s) = - \Big[\frac{R_f}{R_1}X_1(s)+\frac{R_f}{R_2}X_2(s)+\frac{R_f}{R_3}X_3(s) \Big]

Observe que las ganancias del sumador son siempre negativas, hay una inversión de signo en cada señal de entrada.

Y(s) = K_1 X_1(s)+K_2 X_2(s)+K_3 X_3(s)


Ejemplo 1. Implementación con Op-Amp

Referencia: Lathi 4.25 p401

Realizar la implementación de un sistema dado por la función de transferencia H(s)

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

El diagrama de bloques de la función de transferencia H(s) es,

Se agrupan algunos elementos como sumadores y sus factores de multiplicación. Para referencia se etiqueta cada punto como señal W(s) en cada punto donde el orden del exponente de ‘s’ es diferente.

Se considera la inversión de signo de la señal de entrada por la configuración del amplificador operacional y el factor K de cada rama a usar.

Se identifica el tipo de op amp a usar y se establecen los valores de los resistores en múltiplos de 10KΩ y los capacitores en el orden de 10 μF.


Ejemplo 1. Desarrollo analítico

Para revisar el comportamiento del circuito, en caso de implementarlo con OpAmps, el resultado de la función de transferencia para el impulso usando el algoritmo de la sección LTI CT Laplace – Ejercicio resuelto para Y(s)=H(s)X(s) con Sympy-Python

 H(s) = P(s)/Q(s):
   2*s + 5   
-------------
 2           
s  + 4*s + 10
 H(s) en factores:
   2*s + 5   
-------------
 2           
s  + 4*s + 10

H(s) parámetros cuadraticos:
(2*s + 5)/(s**2 + 4*s + 10) : 
 {'A': 2.0, 'B': 5.0, 'a': 2.0, 'c': 10.0,
  'r': 2.041241452319315, 'b': 2.449489742783178,
  'theta': -0.20135792079033082}

 h(t) :
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             

polosceros:
Q_polos : {-2 - sqrt(6)*I: 1, -2 + sqrt(6)*I: 1}
P_ceros : {-5/2: 1}

Estabilidad de H(s):
 n_polos_real : 0
 n_polos_imag : 2
 enRHP : 0
 unicos : 0
 repetidos : 0
 asintota : estable

 X(s): 
1

Respuesta entrada cero ZIR H(s) y condiciones iniciales
term_cero : 0
ZIR :
0
yt_ZIR :
0

 ZSR respuesta estado cero:
ZSR :
   2*s + 5   
-------------
 2           
s  + 4*s + 10

 ZSR_Qs2 :
 (2*s + 5)/(s**2 + 4*s + 10) :
  {'A': 2.0, 'B': 5.0, 'a': 2.0, 'c': 10.0,
   'r': 2.041241452319315, 'b': 2.449489742783178,
   'theta': -0.20135792079033082}
yt_ZSR :
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             

 Y(s)_total = ZIR + ZSR:
   2*s + 5   
-------------
 2           
s  + 4*s + 10

 y(t)_total = ZIR + ZSR:
/  ___  -2*t    /  ___  \                       \             
|\/ 6 *e    *sin\\/ 6 *t/      -2*t    /  ___  \|             
|------------------------ + 2*e    *cos\\/ 6 *t/|*Heaviside(t)
\           6                                   /             
>>>

donde la gráfica de polos muestra que se encuentran todos del lado izquierdo del plano

OpAmpEj01 H(s) Polos

También se muestra la respuesta al impulso h(t) del circuito

Las respuestas fueron obtenidas al usar como bloque de entrada,

# H(s) respuesta impulso
Ps = 2*s+5
Qs = s**2 + 4*s + 10
Hs = Ps/Qs

# X(s) Señal de entrada
xt = d

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [0, 0] # estado cero no se usan

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

4.4.4 LTI Laplace – Ejercicios con circuitos eléctricos y transformadas de Laplace

Ejemplos de solución a circuitos eléctricos que en el planteamiento se usan ecuaciones diferenciales ordinarias y se desarrollan usando la Transformada de Laplace. En el ejercicio 3 se desarrolla una malla que genera un sistema de ecuaciones EDO y se resuelve con transformadas de Laplace.

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejercicio 1. Circuito RLC, fuente DC e interruptor con transformada de Laplace

Referencia: Lathi Ej.4.13 p365

En el circuito de la figura, el interruptor se encuentra cerrado mucho tiempo antes de t=0. Cuando se abre en un instante, encuentre la corriente del inductor y(t) para t≥0.

El interruptor se encuentra cerrado por mucho tiempo, la corriente por el inductor es 2 A y el voltaje del capacitor es 10 V, pues el inductor en DC opera como un conductor sin resistencia y el capacitor se encuentra completamente cargado.

Cuando se abre el interruptor, el circuito se equivalente es el mostrado en la derecha, La corriente inicial del inductor es y(0)=2 y el voltaje inicial del capacitor Vc(0)=10.

El voltaje en la entrada es 10 V, empezando en t=0 y puede ser representado como 10 μ(t) para simplificar la representación de la fuente DC y expresar que antes de t=0 se aplica las condiciones iniciales dadas.mientras se den las condiciones iniciales en t=0 solo será necesario saber que la corriente en t≥0 para determinar la respuesta en t≥0 .

La ecuación del circuito en luego de abrir el interruptor es:

\frac{\delta}{\delta t}y(t) + 2 y8t) + 5 \int_{-\infty}^{t}y(\tau) \delta \tau = 10 \mu(t)

Las condiciones iniciales se aplican como:

\frac{\delta}{\delta t}y(t) = sY(s) - y(0^-)= sY(s)-2 \int_{-\infty}^{t} y(\tau) \delta \tau = \frac{1}{s}Y(s) + \frac{1}{s} \int_{-\infty}^{0^-}y(\tau) \delta \tau

y(t) es la corriente del capacitor, por lo que el integral es:

\int_{-\infty}^{0^-}y(\tau) \delta \tau = q_C(0^-) = CV_c (0^-) = \frac{1}{5}10 = 2

por lo que la parte del capacitor es:

\int_{-\infty}^{t}y(\tau) \delta \tau = \frac{1}{s}Y(s) + \frac{2}{s}

La transformada de Laplace de la ecuación integro diferencial del circuito RLC al usar los resultados se reescribe como:

sY(s)-2+2Y(s)+ 5\frac{1}{s}Y(s) + 5\frac{2}{s} = \frac{10}{s} sY(s)+2Y(s)+ 5\frac{1}{s}Y(s) = \frac{10}{s} +2 - 5\frac{2}{s} \Big[ s+2+\frac{5}{s}\Big]Y(s) =2 \Big[ s^2+2s+5\Big]Y(s) = 2 s Y(s) = \frac{2 s}{s^2+2s+5}

El polinomio del denominador tiene raíces complejas

>>> import  sympy as sym
>>> s = sym.Symbol('s')
>>> Qs = s**2+2*s+5
>>> sym.roots(Qs)
{-1 - 2*I: 1, -1 + 2*I: 1}

por lo que es conveniente usar la forma cuadrática de la transformada de Laplace, donde A= 1, B=0, a=1, c=5, siendo,

r=\sqrt{\frac{20}{4}} = \sqrt{5} b=\sqrt{c-a^2} = 2 \theta=\tan^{-1} \Big(\frac{2}{4} \Big) = 0.46364 rad = 26.56^{o}

usando la tabla de transformadas,

y(t) = \sqrt{5} e^{-t} cos (2t+0.46364) \mu (t)

El gráfico de polos de la función de transferencia en el dominio s:

LTIC Laplace Circuito Eléctrico 01 Polos H(s)

que generan una respuesta de salida y(t)

LTI Laplace Circuito Electrico 01 xh_y

si la entrada x(t) es un impulso unitario, la respuesta es la misma que h(t)

La respuesta del algoritmo obtenida es:

H(s) = P(s)/Q(s):
    2*s     
------------
 2          
s  + 2*s + 5
 H(s) en factores:
    2*s     
------------
 2          
s  + 2*s + 5

H(s) parámetros cuadraticos:
2*s/(s**2 + 2*s + 5) : {'A': 2.0, 'B': 0, 'a': 1.0,
 'c': 5.0, 'r': 2.23606797749979, 'b': 2.0,
 'theta': 0.4636476090008061}

 h(t) :
/   -t               -t         \             
\- e  *sin(2*t) + 2*e  *cos(2*t)/*Heaviside(t)

polosceros:
Q_polos : {-1 - 2*I: 1, -1 + 2*I: 1}
P_ceros : {0: 1}

Estabilidad de H(s):
 n_polos_real : 0
 n_polos_imag : 2
 enRHP : 0
 unicos : 0
 repetidos : 0
 asintota : estable

 X(s): 
1

Respuesta entrada cero ZIR H(s) y condiciones iniciales
term_cero : 0
ZIR :
0
yt_ZIR :
0

 ZSR respuesta estado cero:
ZSR :
    2*s     
------------
 2          
s  + 2*s + 5

 ZSR_Qs2 :
 2*s/(s**2 + 2*s + 5) :
  {'A': 2.0, 'B': 0, 'a': 1.0, 'c': 5.0,
 'r': 2.23606797749979, 'b': 2.0,
 'theta': 0.4636476090008061}
yt_ZSR :
/   -t               -t         \             
\- e  *sin(2*t) + 2*e  *cos(2*t)/*Heaviside(t)

 Y(s)_total = ZIR + ZSR:
    2*s     
------------
 2          
s  + 2*s + 5

 y(t)_total = ZIR + ZSR:
/   -t               -t         \             
\- e  *sin(2*t) + 2*e  *cos(2*t)/*Heaviside(t)
>>>

Instrucciones en Python

Usando los bloques desarrollados en la Unidad 4 Sistemas LTI – Laplace  y las funciones resumidas como telg1001.py que pueden ser usados en cada pregunta.

# Y(s) Respuesta total con entada cero y estado cero
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)
d = sym.DiracDelta(t)
u = sym.Heaviside(t)

# H(s) respuesta impulso
Ps = 2*s
Qs = s**2 + 2*s + 5
Hs = Ps/Qs

# X(s) Señal de entrada
xt = d

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [0, 0] # estado cero no se usan

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_s(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_s_parametros(Hs_fc)

polosceros = fcnm.busca_polosceros(Hs)
Q_polos = polosceros['Q_polos']
P_ceros = polosceros['P_ceros']

estable = fcnm.estabilidad_asintotica_s(Q_polos)

# H(t) respuesta al impulso
ht = 0*s
term_suma = sym.Add.make_args(Hs)
for term_k in term_suma:
    ht_k = sym.inverse_laplace_transform(term_k,s,t)
    # simplifica log(exp()) ej: e**(-2s)/(s**2)
    if ht_k.has(sym.log):
        ht_k = sym.simplify(ht_k,inverse=True)
    ht  = ht + ht_k
lista_escalon = ht.atoms(sym.Heaviside)
ht = sym.expand(ht,t) # terminos suma
ht = sym.collect(ht,lista_escalon)

# PROCEDIMIENTO Respuesta ZIR, ZSR
Xs = fcnm.laplace_transform_suma(xt)

# ZIR_s respuesta entrada cero de s
sol_ZIR = fcnm.respuesta_ZIR_s(Hs,cond_inicio)
ZIR = sol_ZIR['ZIR']
yt_ZIR = sol_ZIR['yt_ZIR']

# ZSR respuesta estado cero, Y(s) a entrada X(s)
sol_ZSR = fcnm.respuesta_ZSR_s(Hs,Xs)
ZSR = sol_ZSR['ZSR']
yt_ZSR = sol_ZSR['yt_ZSR']

# Respuesta total Y(s) y y(t)
Ys = ZIR + ZSR
Ys = fcnm.apart_exp(Ys)
yt = yt_ZIR + yt_ZSR
lista_escalon = yt.atoms(sym.Heaviside)
yt = sym.collect(yt,lista_escalon)

# SALIDA
print(' H(s) = P(s)/Q(s):')
sym.pprint(Hs)
print(' H(s) en factores:')
sym.pprint(Hs_fc)
if len(Hs_Qs2)>0:
    print('\nH(s) parámetros cuadraticos:')
    fcnm.print_resultado_dict(Hs_Qs2)

print('\n h(t) :')
sym.pprint(ht)

print('\npolosceros:')
fcnm.print_resultado_dict(polosceros)

print('\nEstabilidad de H(s):')
for k in estable:
    print('',k,':',estable[k])

print('\n X(s): ')
sym.pprint(Xs)
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')

if not(sol_ZIR == sym.nan): # existe resultado
    fcnm.print_resultado_dict(sol_ZIR)
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

print('\n ZSR respuesta estado cero:')
fcnm.print_resultado_dict(sol_ZSR)

print('\n Y(s)_total = ZIR + ZSR:')
sym.pprint(Ys)
print('\n y(t)_total = ZIR + ZSR:')
sym.pprint(yt)

# Graficas polos, H(s), con polos h(t) --------
figura_s  = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,muestras,f_nombre='H')
figura_ht = fcnm.graficar_ft(ht,t_a,t_b,muestras,f_nombre='h')
# GRAFICAS y(t),x(t),h(t) ---------------------
figura_ft = fcnm.graficar_xh_y(xt,ht,yt,t_a,t_b,muestras)
plt.show()

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejemplo 2. Circuito RLC con transformada de Laplace

Referencia: Lathi Ejemplo 4.17 p375

Realice el análisis de la corriente i(t) en el circuito mostrado en la figura, suponiendo que todas las condicioes iniciales son cero.

El primer paso es representar el circuito en el dominio de la frecuencia (s), mostrado a la derecha de la figura. Se representan los voltajes y corrientes usando la transformada de Laplace.

El voltaje 10μ(t) se representa como 10/s y la corriente i(t) como I(s). Todos los elementos del circuito se muestran con su respectiva impedancia. El inductor de 1 henrio se representa por s, el capacitor de 1/3 faradio se representa por 2/s. El resistor de 3 ohms es solo 3.

El voltaje en el circuito V(s) en cualquier elemento es I(s) por su impedancia Z(s).

La impedancia del circuito para el lazo es:

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

el voltaje de entrada es V(s)=10/s, por lo que la fórmula básica de corriente es:

I(s) =\frac{V(s)}{Z(s)} = \frac{10/s}{(s^2+3s+2)/s} I(s) = \frac{10}{(s^2+3s+2)}

usando las raíces del denominador,

I(s) =\frac{10}{(s+1)(s+2)}

aplicando fracciones parciales:

I(s) = \frac{10}{s+1} -\frac{10}{s+2}

Aplicando la transformada inversa,

i(t) = 10*(e^{-t}-e^{-2t}) \mu (t)

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]


Ejemplo 3. Circuito RC y RL, fuente DC e interruptor con transformada de Laplace

Referencia: Lathi ejemplo 4.18 p378

El interruptor está en posición cerrada por un largo tiempo antes de que sea abierto en t=0. Encuentre las corrientes y1(t) y y2(t) para t>0

Al revisar el circuito se observa que cuando el interruptor esta cerrado y se han alcanzado las condiciones de estado estable, el voltaje del capacitor Vc=16 V y la corriente del inductor y2= 4 .

Las fuentes de 20V y 4V se pueden también ordenar en forma equivalente

Cuando se abre el interruptor en t=0, las condiciones iniciales son Vc(0)=16 y y2(0)= 4 como se muestra en la versión con transformada de Laplace en la figura anterior. La ecuación diferencial del circuito para suma de voltajes en la malla en el lazo 1 y la tabla de propiedades de la transformada de Laplace para la integración aplicada al capacitor,

-\frac{20}{s} + \frac{1}{s}\Big[Y_1(s)+ \int_{-\infty}^{0^-} y_1(\tau) \delta \tau \Big] + \frac{1}{5} \Big[Y_1(s)-Y_2(s) \Big] = 0 -\frac{20}{s} + \frac{1}{s}\Big[Y_1(s) + 16 \Big] + \frac{1}{5} Y_1(s) - \frac{1}{5}Y_2(s) = 0 -\frac{20}{s} + \frac{1}{s}Y_1(s) +\frac{16}{s} + \frac{1}{5}Y_1(s)-\frac{1}{5}Y_2(s) = 0 \Big[ \frac{1}{s}+ \frac{1}{5}\Big] Y_1(s)-\frac{1}{5}Y_2(s) = \frac{4}{s} \frac{1}{5}\Big[ \frac{5+s}{s}\Big] Y_1(s)-\frac{1}{5}Y_2(s) = \frac{4}{s} \Big[ \frac{5+s}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s}
>>> import sympy as sym
>>> s = sym.Symbol('s')
>>> ecuacion1 = -20/s+(1/s)*(Y1+16)+(1/5)*Y1-(1/5)*Y2
>>> ecuacion1.expand()
0.2*Y1 + Y1/s - 0.2*Y2 - 4/s

El voltaje inicial del capacitor se puede representar por un voltaje en serie 16/s y la corriente inicial del inductor de 4 A que representa una fuente de valor VL = Ly2(0) = 1/2(4) = 2. La ecuación diferencial del lazo 2 y usando la transformada de Laplace para la derivada con condiciones iniciales para el inductor,

\frac{1}{5}\Big[Y_2(s)-Y_1(s)\Big] + 1 Y_2(s) + \frac{1}{2} \Big[ sY_2(s) - Y_2(0^-)\big] = 0 \Big[\frac{1}{5}+ 1\Big] Y_2(s) -\frac{1}{5}Y_1(s) + \frac{1}{2} \Big[ sY_2(s) - 4 \Big] = 0 \frac{6}{5} Y_2(s) -\frac{1}{5}Y_1(s) + \frac{1}{2} sY_2(s) - 4\frac{1}{2} = 0 -\frac{1}{5}Y_1(s) +\Big[ \frac{6}{5} + \frac{1}{2} s\Big] Y_2(s) = 2 -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10

el resultado obtenido com Sympy, como una instrucción adicional a la anterior,

>>> ecuacion2 = (1/5)*(Y2-Y1)+Y2+(1/2)*(s*Y2-4)
>>> ecuacion2.expand()
-0.2*Y1 + 0.5*Y2*s + 1.2*Y2 - 2.0

Con lo que hay que resolver el sistema de ecuaciones:

\begin{cases} \Big[ \frac{s+5}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s} \\ -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10 \end{cases}

multipicando la primera ecuación por (s/(5+s)) y sumando con la segunda

\begin{cases} Y_1(s)-\frac{s}{s+5}Y_2(s) = \frac{20}{s}\frac{s}{s+5} \\ -Y_1(s) +\Big[ \frac{12+5s}{2}\Big] Y_2(s) = 10 \end{cases} -\frac{s}{s+5}Y_2(s) + \frac{12+5s}{2}Y_2(s) = \frac{20}{s+5} + 10 \Big[ -\frac{s}{s+5} + \frac{12+5s}{2} \Big] Y_2(s) = 10\Big[\frac{2+(s+5)}{s+5}\Big] \Big[ \frac{-2s+(s+5)(12+5s)}{2(s+5)}\Big] Y_2(s) = 10\frac{2+s+5}{s+5} (-2s+12s+5s^2 +60+25s) Y_2(s) = 20(s+7) (5s^2 +35s +60) Y_2(s) = 20(s+7) 5(s^2 +7s +12) Y_2(s) = 20(s+7) Y_2(s) = 4\frac{s+7}{s^2 +7s +12}

las raices del denominador son s=-4 y s=-3

Y_2(s) = 4\frac{s+7}{(s+3)(s+4)}

usando fracciones parciales y el método de Heaviside,

k_1 = 4\frac{s+7}{\cancel{(s+3)}(s+4)} \Big|_{s=-3} = 4\frac{(-3)+7}{(-3+4)} = 4\frac{4}{1} = 16 k_2 = 4\frac{s+7}{(s+3)\cancel{(s+4)}}\Big|_{s=-4} = 4\frac{(-4)+7}{(-4+3)}=4\frac{3}{-1} = - 12

reescribiendo en fracciones parciales

Y_2(s) = \frac{16}{s+3} -\frac{12}{s+4}

usando la tabla de transformadas de Laplace:

y_2(t) = (16 e^{-3t}-12e^{-4t}) \mu (t)

de forma semejante se puede resolver para Y1(s) al reemplazar el resultado en la primera ecuación del sistema de ecuaciones,

\Big[ \frac{s+5}{s}\Big] Y_1(s)-Y_2(s) = \frac{20}{s} \Big[ \frac{s+5}{s}\Big] Y_1(s)-\Big[\frac{16}{s+3} -\frac{12}{s+4}\Big] = \frac{20}{s} Y_1(s) = \frac{s}{s+5} \Big[\frac{16}{s+3} -\frac{12}{s+4} + \frac{20}{s}\Big] Y_1(s) = \frac{s}{s+5}\frac{16}{s+3} - \frac{s}{s+5}\frac{12}{s+4} + \frac{s}{s+5}\frac{20}{s} Y_1(s) = 16\frac{s}{(s+5)(s+3)} - 12\frac{s}{(s+5)(s+4)} + \frac{20}{s+5}

realizando las fracciones parciales con método de Heaviside para los dos primeros términos,

>>> import sympy as sym
>>> s = sym.Symbol('s')
>>> ya = 16*s/((s+5)*(s+3))
>>> ya.apart()
40/(s + 5) - 24/(s + 3)
>>> yb = -12*s/((s+5)*(s+4))
>>> yb.apart()
-60/(s + 5) + 48/(s + 4)
>>> 
Y_1(s) = \frac{40}{s+5} -\frac{24}{s+3} - \frac{60}{s+5}+\frac{48}{s+4} + \frac{20}{s+5} Y_1(s) = -\frac{24}{s+3} +\frac{48}{s+4}

aplicando desde la tabla la transformada inversa de Laplace,

y_1(t) = (-24e^{-3t} +48e^{-4t} ) \mu (t)

Instrucciones con Python
la solución del sistema de ecuaciones con transformadas de Laplace se puede realizar con Sympy

Y1 (s):
  24*s + 48  
-------------
 2           
s  + 7*s + 12

Y1 (s) en fracciones parciales:
  48      24 
----- - -----
s + 4   s + 3
Y2 (s):
   4*s + 28  
-------------
 2           
s  + 7*s + 12

Y2 (s) en fracciones parciales:
    12      16 
- ----- + -----
  s + 4   s + 3
>>> 
# Sistemas de ecuaciones con Sympy
import sympy as sym

# INGRESO
s = sym.Symbol('s')
[Y1, Y2] = sym.symbols('Y1 Y2')

ecuacion1 = ((s+5)/s)*Y1 - Y2 -20/s
ecuacion2 = -Y1 + +((12+5*s)/2)*Y2 -10

variables = [Y1,Y2]

# PROCEDIMIENTO
respuesta = sym.solve([ecuacion1,
                       ecuacion2],
                       variables)

# SALIDA
for cadarespuesta in respuesta:
    print(str(cadarespuesta) +'(s):')
    sym.pprint(respuesta[cadarespuesta])
    print('')
    print(str(cadarespuesta)+'(s) en fracciones parciales:') 
    sym.pprint(respuesta[cadarespuesta].apart())

[Ej 1. RLC fuente DC e interruptor] [ Ej 2. RLC con Laplace ] [Ej 3. RC, RL, interruptor ]

4.4.3 LTI Laplace – Ejercicios resueltos de sistemas Lineales contínuos LTI

Ejemplos, ejercicios resueltos con Transformada de Laplace aplicados en sistemas Lineales contínuos LTI.

Ejercicios 1ra Evaluación


1Eva2012TI_T1 LTI CT Polos y ceros de H(s) con respuesta DC conocida

4.4.2 LTI Laplace – Ejercicio para Y(s)=H(s)X(s) con constantes símbolo con Sympy-Python

Los ejercicios incluyen constantes tipo símbolo, que por simplicidad de la ecuación en la resolución se usan solo las operaciones básicas para la transformada Inversa de Laplace.

Ejemplo 1.  Y(s) dado H(s) y X(s) con una constante α

Referencia: Oppenheim ejemplo 9.37 p718 pdf749, semejante al ejercicio sistema-modelo LTI Lathi ejemplo 1.16 p111. Oppenheim problema 2.61c p164 pdf195

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

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

que en el dominio ‘s’ se escribe como

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

Que es semejante a lo realizado en el ejemplo 2 de H(s) Diagramas de bloques, dominio ‘s’ y ecuaciones diferenciales, de donde se tiene que:

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

Suponga que la señal de entrada es x(t) = α μ(t) , usando la tabla de transformadas de Laplace,  la entrada del sistema es,

X(s) = \frac{\alpha}{s} = \frac{\alpha}{s+0}

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+1)(s+2)} \Big]\Big[\frac{\alpha}{s+0}\Big]

usando fracciones parciales con Python para encontrar el equivalente, la expresión se puede separar en,

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

Luego, aplicando la transformada inversa de Laplace de 1/(s+a) que es una exponencial e-at μ(t), se tiene la respuesta en el dominio del tiempo.

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

Instrucciones con Sympy-Python

Se crea la expresión H(s) y X(s) en forma simbólica usando las variables, s y a.

Dado que la expresión es tipo polinomio, y que no existen exponenciales, el desarrollo del algoritmo puede realizarse con operaciones básicas realizadas en las  secciones anteriores. Añadiendo la operación para encontrar la salida Y(s)

# Señal de salida Y(s)
Ys = Hs*Xs
Ys_parcial = Ys.apart(s)

Además de realizar la transformada inversa de Laplace en cada término de Y(s).  Con lo que se obtiene el resultado:

Ys:
       a        
----------------
  / 2          \
s*\s  + 3*s + 2/

 Ys = Hs*Xs :
    a         a      a 
--------- - ----- + ---
2*(s + 2)   s + 1   2*s

 y(t) :
                                         -2*t             
a*Heaviside(t)      -t                a*e    *Heaviside(t)
-------------- - a*e  *Heaviside(t) + --------------------
      2                                        2          
>>> sym.pprint(yt.simplify())
  / 2*t      t    \  -2*t             
a*\e    - 2*e  + 1/*e    *Heaviside(t)
--------------------------------------
                  2                   
>>>   

Instrucciones con Sympy-Python

# Y(s) = H(s)*X(s)
# Oppenheim ejemplo 9.37 p718/pdf746
import sympy as sym

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)
a = sym.Symbol('a', real=True)

Hs = 1/(s**2+3*s+2)
Xs = a/s

# PROCEDIMIENTO
# Señal de salida Y(s)
Ys = Hs*Xs
Ysp = Ys.apart(s)

yt = sym.inverse_laplace_transform(Ysp,s,t)

# SALIDA
print('Ys:')
sym.pprint(Ys)
print('\n Ys = Hs*Xs :')
sym.pprint(Ysp)
print('\n y(t) :')
sym.pprint(yt)

Ejemplo 2. Y(s) con condiciones iniciales tipo constante simbolo

Referencia: Oppenheim ejemplo 9.38 p719

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

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

sea x(t) = αμ(t)

X(s) = \frac{\alpha}{s} = \frac{\alpha}{s+0}

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

\frac{\delta^2}{\delta t^2}y(t) \Rightarrow s^2Y(s) - y(0^{-})s - y'(0^{-}) = s^2Y(s) - \beta s - \gamma \frac{\delta}{\delta t}y(t) \Rightarrow sY(s) - y(0^{-}) = sY(s) - \beta

resultados que se insertan en la ecuación original de derivadas de t:

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

que se reordenan para despejar Y(s),

[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)} Y(s) = \frac{\beta (s+3)}{(s+1)(s+2)} + \frac{\gamma}{(s+1)(s+2)} + \frac{\alpha}{s(s+1)(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

Instrucciones con Sympy-Python

En este caso, se ingresan las propiedades de diferenciación para Laplace, aunque pueden crearse a partir de una fórmula mas general en la siguiente sección.

La ecuación se define como una expresión de (lado_izquierdo)=(lado_derecho) contruida a partir la ecuación del problema incluyendo las condiciones iniciales diferentes de cero.

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 sustituir los valores de a,b y g, se usa un diccionario {} con las parejas de {variable: valor}.

ecuacion:
   2                                 a
Y*s  + 3*Y*s + 2*Y - b*s - 3*b - g = -
                                     s

 Y(s): 
       2              
a + b*s  + 3*b*s + g*s
----------------------
     / 2          \   
   s*\s  + 3*s + 2/   

 Y(s) en Fracciones parciales
 a    a - 2*b - 2*g   a - 2*b - g
--- + ------------- - -----------
2*s     2*(s + 2)        s + 1   

sustituye valores de a, b y g
 Ys:
  3       1     1
----- - ----- + -
s + 2   s + 1   s

 y(t) con constantes simbolo 
/   2*t                                     t\  -2*t             
\a*e    + a - 2*b - 2*g + 2*(-a + 2*b + g)*e /*e    *Heaviside(t)
-----------------------------------------------------------------
                                2                                

sustituye valores de a, b y g

 y(t): 
/ 2*t    t    \  -2*t             
\e    - e  + 3/*e    *Heaviside(t)

 y(t) expande expresion
                -t                   -2*t             
Heaviside(t) - e  *Heaviside(t) + 3*e    *Heaviside(t)
>>> 

Compare los resultados con lo expuesto en la parte teórica.

# Y(s) con condiciones iniciales como constantes
# Oppenheim ejemplo 9.38 p719
import numpy as np
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
s = sym.Symbol('s')
a,b,g = sym.symbols('a b g')
Y = sym.Symbol('Y')

# propiedad de diferenciacion
d2Y = (s**2)*Y - b*s - g
dY  = s*Y - b

Xs = a/s    # señal de entrada

# PROCEDIMIENTO
# ecuacion del sistema
Lizq = d2Y + 3*dY + 2*Y # Lado izquierdo
Lder = Xs
ecuacion = sym.Eq(Lizq, Lder)

# despeja la ecuacion para Y(s)
Ys = sym.solve(ecuacion,Y)[0]

Yspc = sym.apart(Ys,s) # incluye constante simbolo

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

# y(t) con inversa transformada Laplace
ytc = sym.inverse_laplace_transform(Yspc,s,t)
yt = sym.inverse_laplace_transform(Ysp,s,t)

# SALIDA
print('ecuacion:')
sym.pprint(ecuacion)
print('\n Y(s): ')
sym.pprint(Ys)
print('\n Y(s) en Fracciones parciales')
sym.pprint(Yspc)
print('\nsustituye valores de a, b y g')
print(' Ys:')
sym.pprint(Ysp)
print('\n y(t) con constantes simbolo ')
sym.pprint(ytc)
print('\nsustituye valores de a, b y g')
print('\n y(t): ')
sym.pprint(yt)
print('\n y(t) expande expresion')
sym.pprint(sym.expand(yt))

4.4.1 LTI Laplace – Ejercicio resuelto para Y(s)=H(s)X(s) con Sympy-Python

Se presenta un ejercicio que integra los temas tratados en la Unidad 4 para el análisis del Sistema con Transformada de Laplace. El algoritmo se muestra como una referencia para adaptar en los ejercicios integradores usando las funciones resumidas en telg1001.py .

Referencia: Lathi Ejemplo 4.12 p361

Resuelva la ecuación diferencial de segundo orden

(D^2 + 5D + 6) y(t) = (D+1) x(t)

para las condiciones iniciales y(0-)=2 y y'(0-)=1, cuando se presenta la entrada

x(t) = e^{-4t} \mu (t)

Desarrollo

La ecuación diferencial del sistema, interpretando los operadores D por derivadas se muestra a continuación.

\frac{\delta^2}{\delta t^2}y(t) + 5\frac{\delta}{\delta t}y(t) + 6y(t) = \frac{\delta^2}{\delta t^2}x(t) + x(t)

La ecuación se representa en el dominio s de las transformadas de Laplace para usarlas en el algoritmo:

s^2 Y(s) + 5s Y(s) + 6 Y(s) = s X(s) + X(s) (s^2 + 5s + 6)Y(s) = (s+1) X(s)

con lo que se escribe la función de transferencia H(s) que se ingresa para el análisis del algoritmo.

H(s) = \frac{X(s)}{Y(s)} = \frac{s+1}{s^2 + 5s + 6}

Para el caso de x(t) se usa la transformada de Laplace con el algoritmo para obtener X(s), mostrando que se puede realizar directamente para términos simples. Si la expresión es mas elaborada se usaría las funciones realizadas para f(t) con Sympy-Python

Con el algoritmo se obtienen los siguientes resultados, con lo que se puede analizar la estabilidad del sistema usando la gráfica de polos. No se observan polos en el lado derecho del plano RHP, por lo que el sistema es asintóticamente estable, BIBO-estable.

LTIC Laplace Algoritmo Integrador 01 polos

la gráfica de H(s) junto a los polos:

LTIC Laplace Algoritmo Integrador 01 polos Hs

la gráfica de h(t) como respuesta al impulso,

LTIC Laplace Algoritmo Integrador 01 ht

la gráfica de la señal de entrada x(t) y la respuesta y(t) es:

LTIC Laplace Algoritmo Integrador 01 xh_y

Las expresiones para Y(s), H(s), análisis de estabilidad que entrega el algoritmo son:

 H(s) = P(s)/Q(s):
  2       1  
----- - -----
s + 3   s + 2
 H(s) en factores:
     s + 1     
---------------
(s + 2)*(s + 3)

 h(t) :
/   -2*t      -3*t\             
\- e     + 2*e    /*Heaviside(t)

polosceros:
Q_polos : {-2: 1, -3: 1}
P_ceros : {-1: 1}

Estabilidad de H(s):
 n_polos_real : 2
 n_polos_imag : 0
 enRHP : 0
 unicos : 0
 repetidos : 0
 asintota : estable

 X(s): 
  1  
-----
s + 4

Respuesta entrada cero ZIR H(s) y condiciones iniciales
term_cero : 2*s + 11
ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2
yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)

 ZSR respuesta estado cero:
ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)
yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             

 Y(s)_total = ZIR + ZSR:
      3         3         13   
- --------- - ----- + ---------
  2*(s + 4)   s + 3   2*(s + 2)

 y(t)_total = ZIR + ZSR:
/    -2*t                -4*t\             
|13*e          -3*t   3*e    |             
|-------- - 3*e     - -------|*Heaviside(t)
\   2                    2   /             
>>>

Con los datos de H(s) y H(s) con fracciones parciales, se presentan dos formas de realizar el diagrama de bloques. El primer diagrama corresponde a H(s), mientras qu el segundo corresponde a la interpretación de las fracciones parciales.

Instrucciones Python

Usando los bloques desarrollados en la Unidad 4 Sistemas LTI – Laplace  y las funciones resumidas como telg1001.py que pueden ser usados en cada pregunta, se agrupan y tiene como resultado:

# Y(s) Respuesta total con entada cero y estado cero
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)
d = sym.DiracDelta(t)
u = sym.Heaviside(t)

# H(s) respuesta impulso
Ps = s+1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 1+0*s cuando es constante

# X(s) Señal de entrada
xt = sym.exp(-4*t)*u

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [1, 2] # estado cero no se usan

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_s(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_s_parametros(Hs_fc)

polosceros = fcnm.busca_polosceros(Hs)
Q_polos = polosceros['Q_polos']
P_ceros = polosceros['P_ceros']

estable = fcnm.estabilidad_asintotica_s(Q_polos)

# H(t) respuesta al impulso
ht = 0*s
term_suma = sym.Add.make_args(Hs)
for term_k in term_suma:
    ht_k = sym.inverse_laplace_transform(term_k,s,t)
    # simplifica log(exp()) ej: e**(-2s)/(s**2)
    if ht_k.has(sym.log):
        ht_k = sym.simplify(ht_k,inverse=True)
    ht  = ht + ht_k
lista_escalon = ht.atoms(sym.Heaviside)
ht = sym.expand(ht,t) # terminos suma
ht = sym.collect(ht,lista_escalon)

# PROCEDIMIENTO Respuesta ZIR, ZSR
Xs = fcnm.laplace_transform_suma(xt)

# ZIR_s respuesta entrada cero de s
sol_ZIR = fcnm.respuesta_ZIR_s(Hs,cond_inicio)
ZIR = sol_ZIR['ZIR']
yt_ZIR = sol_ZIR['yt_ZIR']

# ZSR respuesta estado cero, Y(s) a entrada X(s)
sol_ZSR = fcnm.respuesta_ZSR_s(Hs,Xs)
ZSR = sol_ZSR['ZSR']
yt_ZSR = sol_ZSR['yt_ZSR']

# Respuesta total Y(s) y y(t)
Ys = ZIR + ZSR
Ys = fcnm.apart_s(Ys)
yt = yt_ZIR + yt_ZSR
lista_escalon = yt.atoms(sym.Heaviside)
yt = sym.collect(yt,lista_escalon)

# SALIDA
print(' H(s) = P(s)/Q(s):')
sym.pprint(Hs)
print(' H(s) en factores:')
sym.pprint(Hs_fc)
if len(Hs_Qs2)>0:
    print('\nH(s) parámetros cuadraticos:')
    fcnm.print_resultado_dict(Hs_Qs2)

print('\n h(t) :')
sym.pprint(ht)

print('\npolosceros:')
fcnm.print_resultado_dict(polosceros)

print('\nEstabilidad de H(s):')
for k in estable:
    print('',k,':',estable[k])

print('\n X(s): ')
sym.pprint(Xs)
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')

if not(sol_ZIR == sym.nan): # existe resultado
    fcnm.print_resultado_dict(sol_ZIR)
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

print('\n ZSR respuesta estado cero:')
fcnm.print_resultado_dict(sol_ZSR)

print('\n Y(s)_total = ZIR + ZSR:')
sym.pprint(Ys)
print('\n y(t)_total = ZIR + ZSR:')
sym.pprint(yt)

# Graficas polos, H(s), con polos h(t) --------
figura_s  = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,muestras,f_nombre='H')
figura_ht = fcnm.graficar_ft(ht,t_a,t_b,muestras,f_nombre='h')
# GRAFICAS y(t),x(t),h(t) ---------------------
figura_ft = fcnm.graficar_xh_y(xt,ht,yt,t_a,t_b,muestras)
plt.show()

4.4 LTI Laplace – Y(s)=ZIR+ZSR Respuesta entrada cero y estado cero con Sympy-Python

Referencia: Lathi Ejemplo 4.12 p361

Resuelva la ecuación diferencial de segundo orden usando respuesta a entrada cero y respuesta a estado cero.

(D^2 + 5D + 6) y(t) = (D+1) x(t)

para las condiciones iniciales y(0-)=2 y y'(0-)=1, ante la entrada

x(t) = e^{-4t} \mu (t)

El ejercicio muestra nuevamente el concepto sobre la respuesta total del sistema tiene dos componentes:

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]
..


1. Respuesta a entrada cero y condiciones iniciales

Si el sistema bajo observación se no tiene señal de entrada X(s) = 0, la parte de estado cero no aporta componente de salida.

Respuesta
total
= respuesta a
entrada cero ZIR
+ respuesta a
estado cero ZSR

Sin embargo, si existen valores iniciales en los componentes del sistema (circuito con capacitores), se producirá una respuesta a entrada cero que básicamente responde a las condiciones iniciales al sistema.

Referencia: Lathi Ejemplo 4.12 p361

Un sistema LTI tiene la ecuación diferencial de segundo orden expresada con operadores D:

(D^2 + 5D + 6) y(t) = (D+1) x(t)

El sistema tiene las condiciones iniciales y(0-)=2 y y'(0-)=1.

1.1 Desarrollo Analítico para Respuesta a entrada cero y condiciones iniciales

La ecuación se escribe usando la transformada de Laplace:

(s^2 + 5s + 6) Y(s) = (s+1) X(s)

El diagrama básico de la expresión LTI, sin simplificar, se muestra en la figura. La sección de color naranja representa Q(s)  y la sección en verde representa P(s)

Para el análisis de entrada cero X(s)=0, la ecuación se simplifica a

(s^2 + 5s + 6) Y(s) = 0

Las condiciones iniciales se aplican usando la tabla de propiedades de Laplace, donde se muestra que los valores iniciales se añaden a los componentes de primera y segunda derivada:

\frac{\delta}{\delta t} y(t) \Longleftrightarrow sY(s) - y(0^{-}) = sY(s) -2 \frac{\delta^2}{\delta t^2} y(t) \Longleftrightarrow s^2 Y(s) - sy(0^{-}) - y'(0^{-}) = s^2 Y(s) -2s - 1

Se recuerda que para la entrada x(t) = 0, no se crean términos en el lado derecho de la expresión

[s^2 Y(s) -2s - 1] + 5[sY(s) -2] + 6Y(s) = 0

Agrupando los terminos Y(s), se tienen los componentes que no dependen de Y(s) que son los resultantes de las condiciones iniciales.

[s^2 + 5s + 6]Y(s) - (2s+11) = 0 [s^2 + 5s + 6]Y(s) = (2s+11) Y(s)= \frac{2s+11}{s^2 + 5s + 6}

El resultado Y(s) corresponde al término de entrada cero con condiciones iniciales aplicadas.

Aplicando un algoritmo en Python para añadir los componentes de condiciones iniciales, el resultado es:

Respuesta entrada cero ZIR H(s) y condiciones iniciales

 term_cero :
2*s + 11

 ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2

 yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)
>>>

1.2 Desarrollo con Sympy-Python

Se inicia con los resultados de algoritmos anteriores, donde se ingresan las expresiones de Q(s) y P(s) correspondientes a H(s). Las condiciones iniciales  se ingresan en un arreglo en orden de derivada descendente, por ejemplo para  el ejercicio [y'(0), y(0)] = [1., 2.]

Par generar las expresiones de condicion inicial, se genera el término de condición cero o term_cero, compuesto de los términos de orden tomados del vector inicial.

# Y_ZIR(s) Respuesta entrada cero con Transformada de Laplace
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

# H(s) respuesta impulso
Ps = s + 1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 6*sym.exp(-2*s)/s**2

# X(s) Señal de entrada
Xs = 0*s

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [1, 2]

# PROCEDIMIENTO
def respuesta_ZIR_s(Hs,cond_inicio):
    '''respuesta a entrada cero ZIR con H(s) y condiciones de inicio
       Si ZIR es sym.nan existen insuficientes condiciones
    '''
    polosceros = fcnm.busca_polosceros(Hs)
    Q_polos = polosceros['Q_polos']
    
    # coeficientes y grado de Q
    Q = 1+0*s
    for raiz in Q_polos:
        veces = Q_polos[raiz]
        Q = Q*((s-raiz)**veces)
    Q = sym.poly(Q,s)
    
    Q_coef  = Q.all_coeffs() # coeficientes Q
    Q_grado = Q.degree()  # grado polinomio Q
    
    term_cero = 0*s 
    if len(cond_inicio) == Q_grado:
        for j in range(0,Q_grado,1):
            term_orden = 0
            for i in range(j,Q_grado,1):
                term_orden = term_orden + cond_inicio[i]*(s**(i-j))
            term_cero = term_cero + term_orden*Q_coef[j]
        ZIR = term_cero/Q
        ZIR = fcnm.apart_s(ZIR)
        ZIR_Qs2 = fcnm.Q_cuad_s_parametros(ZIR)
    else:
        ZIR = sym.nan # insuficientes condiciones iniciales

    if not(ZIR==sym.nan):
        # entrada_cero en t
        yt_ZIR = sym.inverse_laplace_transform(ZIR,s,t)
        # simplifica log(exp()) ej: e**(-2t)/(s**2)
        if yt_ZIR.has(sym.log):
            yt_ZIR = sym.simplify(yt_ZIR,inverse=True)
        lista_escalon = yt_ZIR.atoms(sym.Heaviside)
        yt_ZIR = sym.expand(yt_ZIR,t) # terminos suma
        yt_ZIR = sym.collect(yt_ZIR,lista_escalon)
        
        sol_ZIR ={'term_cero' : term_cero}
        if len(ZIR_Qs2)>0: # añade si Q es cuadratico
            sol_ZIR['ZIR_Qs2'] = ZIR_Qs2
        sol_ZIR = sol_ZIR | {'ZIR'   : ZIR,
                             'yt_ZIR': yt_ZIR}
    else:
        sol_ZIR = sym.nan
    
    return(sol_ZIR)

# respuesta a entrada_cero ZIR en s
sol_ZIR = respuesta_ZIR_s(Hs,cond_inicio)

# SALIDA - polinomio
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')
if not(sol_ZIR == sym.nan): # existe resultado
    for k in sol_ZIR:
        print('\n',k,':')
        sym.pprint(sol_ZIR[k])
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

El siguiente caso sería realizar el mismo sistema LTI CT, con respuesta Y(s) ante una extrada X(s), conocido también como estado cero.

Para mantenerlo simple, se consideraría condiciones iniciales en cero.

Finalmente se sumarían la entrada cero y estado cero para obtener la respuesta total.

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

..


2. Respuesta a estado cero con x(t) de Sistema con ecuación diferencial de 2do orden

El ejercicio continua  el desarrollo del concepto sobre la respuesta total del sistema tiene dos componentes:

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero
ZSR = h(t) ⊗ x(t)
ZSR(s) = H(S) X(s)

Referencia: Lathi Ejemplo 4.12 p361

Un sistema LTI tiene la ecuación diferencial de segundo orden expresada con operadores D:

(D^2 + 5D + 6) y(t) = (D+1) x(t)

(viene de la sección anterior…) si se somete a la entrada,

x(t) = e^{-4t} \mu (t)

2.1. Desarrollo Analítico para estado cero

Ahora, todos los componentes tienen condición inicial cero y una señal de entrada x(t). El diagrama de bloques del sistema se obtuvo en la sección anterior de «entrada cero»:

En este caso se considera que ningún componente del sistema tiene valor inicial [0., 0.], interpretado por ejemplo como que todos los capacitores estan descargados, la salida corresponde solo a la parte conformada por x(t) y la respuesta al impulso h(t).

para la entrada x(t),

x(t) = e^{-4t} \mu(t) \Longleftrightarrow X(s) = \frac{1}{s+4}

Usando transformadas de Laplace, la ecuación del sistema se convierte en:

[s^2 Y(s)] + 5[sY(s)] + 6Y(s) = (s+1)X(s)

que al aplicar la parte de X(s),

[s^2 Y(s)] + 5[sY(s)] + 6Y(s) = (s+1)\frac{1}{s+4} (s^2 + 5s + 6) Y(s) = \frac{s+1}{s+4} Y(s) = \frac{s+1}{(s+4)(s^2 + 5s + 6)}

que es equivalente a y(t) = h(t)⊗x(t) como convolución en dominio de tiempo o Y(s)=H(s)X(s) en el dominio de s.

Para simplificar la respuesta de la inversa de transformada de Laplace, se separa en fracciones parciales. Al usar el algoritmo en Python se obtiene:

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

Tomando la referencia de la tabla de transformadas de Laplace, fila 5, se tiene:

y(t) = -\frac{1}{2} e^{-2t} \mu(t) + 2 e^{-3t} \mu(t) -\frac{3}{2} e^{-4t} \mu (t)

2.2. Instrucciones con Sympy-Python

Para este caso se añade la transformada de Laplace X(s) de x(t). La expresión μ(t) se representa como Heaviside en Sympy.

Las instrucciones son semejantes a las usadas en algoritmos anteriores, tal como Hs = Ps/Qs, considerando que respuesta a estado cero es ZSR = Hs*Xs

Se refleja el resultado obtenido en el desarrollo analítico,

 ZSR respuesta estado cero:

 ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)

 yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             
>>> 

las instrucciones usadas son:

# Y_ZSR(s) Respuesta estado cero con Transformada Laplace
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

# H(s) respuesta impulso
Ps = s + 1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs

# X(s) Señal de entrada
Xs = 1/(s+4)

# PROCEDIMIENTO
def respuesta_ZSR_s(Hs,Xs):
    '''respuesta a estado cero ZSR con H(s) y X(s)
    '''
    ZSR = Hs*Xs
    ZSR = fcnm.apart_s(ZSR)
    ZSR_Qs2 = fcnm.Q_cuad_s_parametros(ZSR)

    # ZSR respuesta estado cero(t)
    yt_ZSR = sym.inverse_laplace_transform(ZSR,s,t)

    # simplifica log(exp()) ej: e**(-2t)/(s**2)
    if yt_ZSR.has(sym.log):
        yt_ZSR = sym.simplify(yt_ZSR,inverse=True)
    lista_escalon = yt_ZSR.atoms(sym.Heaviside)
    yt_ZSR = sym.expand(yt_ZSR,t) # terminos suma
    yt_ZSR = sym.collect(yt_ZSR,lista_escalon)

    sol_ZSR = {'ZSR' : ZSR}
    if len(ZSR_Qs2)>0:
        sol_ZSR['ZSR_Qs2'] = ZSR_Qs2
    sol_ZSR = sol_ZSR | {'yt_ZSR' : yt_ZSR}
    return(sol_ZSR)

sol_ZSR = respuesta_ZSR_s(Hs,Xs)

# SALIDA - polinomio
print('\n ZSR respuesta estado cero:')
for k in sol_ZSR:
    print('\n',k,':')
    sym.pprint(sol_ZSR[k])

El siguiente bloque consiste en hacer el ejercicio completo del sistema LTIC, considerando condiciones iniciales y una entrada x(t), que da la respuesta total del sistema.

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

..


3. Y(s) respuesta total con condiciones iniciales y entrada x(t)

Finalmente se obtiene el resultado sumando los dos componentes desarrollados en  los numerales anteriores

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero

3.1 Desarrollo Analítico

Se incluye para mostrar que la respuesta total Y(s) = entrada_cero+estado_cero

Para resolver con transformadas de Laplace, la ecuación se expresa como :

\frac{\delta ^2}{\delta t^2} y(t) + 5 \frac{\delta}{\delta t} y(t) + 6 y(t) = \frac{\delta}{\delta t}x(t) +x(t)

haciendo y(t) ⇔Y(s) y usando la tabla de propiedades

\frac{\delta}{\delta t} y(t) \Longleftrightarrow sY(s) - y(0^{-}) = sY(s) -2 \frac{\delta^2}{\delta t^2} y(t) \Longleftrightarrow s^2 Y(s) - sy(0^{-}) - y'(0^{-}) = s^2 Y(s) -2s - 1

para la entrada x(t),

x(t) = e^{-4t} \mu(t) \Longleftrightarrow X(s) = \frac{1}{s+4} \frac{\delta}{\delta t}x(t) \Longleftrightarrow sX(s) - x(0^{-}) = s\frac{1}{s+4} -0 = \frac{s}{s+4}

con lo que la ecuación diferencial en transformada de Laplace se escribe,

[s^2 Y(s) -2s - 1] + 5[sY(s) -2] + 6Y(s) = \frac{s}{s+4}+\frac{1}{s+4}

Agrupando los terminos Y(s) y los términos de la derecha,

[s^2 + 5s + 6]Y(s) - (2s+11) = \frac{s+1}{s+4} [s^2 + 5s + 6]Y(s) = (2s+11) + \frac{s+1}{s+4} Y(s) = \frac{2s+11}{s^2 + 5s + 6} + \frac{s+1}{(s+4)(s^2 + 5s + 6)}

despejando Y(s) se puede observar que la expresión es la suma de,

Y(s) = \text{entrada cero} + \text{estado cero} Y(s) = \frac{(s+4)(2s+11)+s+1}{(s+4)(s^2 + 5s + 6)} Y(s)= \frac{2s^2+20s+45}{(s^2 + 5s + 6)(s+4)} Y(s)= \frac{2s^2+20s+45}{(s+2)(s+3)(s+4)}

expandiendo en fracciones parciales:

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

Finalmente, usando la transformada inversa de Laplace obtener la expresión en el dominio de t,

Y(t)= \Big[ \frac{13}{2}e^{-2t} - 3 e^{-3t} - \frac{3}{2}e^{-4t} \Big] \mu(t)

El ejercicio muestra la facilidad de resolver el sistema usando transformadas de Laplace.

Respuesta entrada cero ZIR H(s) y condiciones iniciales

 term_cero :
2*s + 11

 ZIR :
    5       7  
- ----- + -----
  s + 3   s + 2

 yt_ZIR :
/   -2*t      -3*t\             
\7*e     - 5*e    /*Heaviside(t)

 ZSR respuesta estado cero:

 ZSR :
      3         2         1    
- --------- + ----- - ---------
  2*(s + 4)   s + 3   2*(s + 2)

 yt_ZSR :
/   -2*t                -4*t\             
|  e          -3*t   3*e    |             
|- ----- + 2*e     - -------|*Heaviside(t)
\    2                  2   /             

 Y0(s) total = ZIR + ZSR:
      3         3         13   
- --------- - ----- + ---------
  2*(s + 4)   s + 3   2*(s + 2)

 y(t) total = ZIR + ZSR:
/    -2*t                -4*t\             
|13*e          -3*t   3*e    |             
|-------- - 3*e     - -------|*Heaviside(t)
\   2                    2   /             
>>> 

3.2 Instrucciones Sympy-Python

Las instrucciones para respuesta a entrada cero ZIR y respuesta estado cero ZSR se integran en un algoritmo.

# Y(s) Respuesta total con entrada cero y estado cero
# Qs Y(s) = Ps X(s) ; H(s)=Ps/Qs
# Ejemplo Lathi 4.12 p361
# http://blog.espol.edu.ec/telg1001/
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

# H(s) respuesta impulso
Ps = s + 1
Qs = s**2 + 5*s + 6
Hs = Ps/Qs
#Hs = 6*sym.exp(-2*s)/s**2

# X(s) Señal de entrada
Xs = 1/(s+4)

# condiciones iniciales, [y'(0),y(0)] orden descendente
t0 = 0
cond_inicio = [1, 2]

# PROCEDIMIENTO
# ZIR_s respuesta entrada cero de s
sol_ZIR = fcnm.respuesta_ZIR_s(Hs,cond_inicio)
ZIR = sol_ZIR['ZIR']
yt_ZIR = sol_ZIR['yt_ZIR']

# ZSR respuesta estado cero, Y(s) a entrada X(s)
sol_ZSR = fcnm.respuesta_ZSR_s(Hs,Xs)
ZSR = sol_ZSR['ZSR']
yt_ZSR = sol_ZSR['yt_ZSR']

# Respuesta total Y(s) y y(t)
Ys = ZIR + ZSR
Ys = fcnm.apart_s(Ys)
yt = yt_ZIR + yt_ZSR
lista_escalon = yt.atoms(sym.Heaviside)
yt = sym.collect(yt,lista_escalon)

# SALIDA - polinomio
print('\nRespuesta entrada cero ZIR H(s) y condiciones iniciales')
if not(sol_ZIR == sym.nan): # existe resultado
    for k in sol_ZIR:
        print('\n',k,':')
        sym.pprint(sol_ZIR[k])
else:
    print(' insuficientes condiciones iniciales')
    print(' revisar los valores de cond_inicio[]')

print('\n ZSR respuesta estado cero:')
for k in sol_ZSR:
    print('\n',k,':')
    sym.pprint(sol_ZSR[k])

print('\n Y0(s) total = ZIR + ZSR:')
sym.pprint(Ys)
print('\n y(t) total = ZIR + ZSR:')
sym.pprint(yt)

Respuesta del sistema a: [ entrada cero ZIR ]  [ estado cero ZSR ]  [total ]

4.3.3 H(s) – polos y Estabilidad del sistema con Sympy-Python

Referencia: Lathi 4.3-3 p371, Lathi 2.4-2 p198

En un sistema acotado (bounded) estable, sometido a una señal de entrada en la que se conoce que su salida siempre es la misma, se lo considera un sistema BIBO-estable. (Bounded input, Boundes output)

La función de transferencia H(s) de un sistema es una descripción externa con la forma:

H(s) = \frac{P(s)}{Q(s)}

El polinomio Q(s) representa una descripción interna, en la que si todos los polos de H(s) se encuentran en la parte izquierda del plano (left half plane, LHP), todos los términos h(t) son exponenciales decrecientes y h(t) es absolutamente integrable. En consecuencia el sistema es BIBO-estable.

Se asume que H(s) es una función en la que M ≤ N, siendo M el grado de P(s) y N e grado de Q(s), caso contrario el sistema es es BIBO-inestable.

Laplace estabilidad de un polo real animación

1. Un sistema LTI CT es asintóticamente estable si y solo si todos los polos de su función de transferencia H(s) se encuentran en el lado izquierdo del plano, pudiendo ser simples o repetidos.

2. Un sistema LTI CT es inestable si y solo si se dan una o ambas de las condiciones:

2.1 Existe al menos un polo de H(s) en el lado derecho del plano RHP.
2.2 Existen polos repetidos de H(s) en el eje imaginario.

3. Un sistema LTI CT es marginalmente estable si y solo si, no hay polos de H(s) en el lado derecho del plano RHP y algunos polos NO repetidos sobre en el eje imaginario (parte Real = 0).

Laplace Establilidad Dos Polos Imaginarios

 


polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]
..


Ejemplo 1. H(s) con polos no repetidos, RHP=0

Referencia: Lathi 2.14a p201

Revise la estabilidad BIBO de los siguientes sistemas LTI¨CT descritos por las ecuaciones:

(D+1)(D^2 + 4D +8) y(t) = (D-3) x(t)

Ejemplo 1 Desarrollo analítico

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+1)(s^2 + 4s +8) Y(s) = (s-3) X(s) \frac{Y(s)}{X(s)} = \frac{s-3}{(s+1)(s^2 + 4s +8)}

El grado del polinomio Q es N=3 y es mayor que el grado del polinomio P es M=1.

La expresión de entrada se escribe en el algoritmo como:

Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

Los polos se obtienen como las raíces del denominador Q, necesario para plantear las fracciones parciales y simplificar la expresión como la suma de componentes más simples.

polos: {-1: 1, -2 - 2*I: 1, -2 + 2*I: 1}

Las fracciones parciales de H(s) usando el algoritmo de la sección anterior se obtiene

H(s) = \frac{1}{5}\frac{4}{s+1}+\frac{1}{5} \frac{4s+17}{s^2 + 4s +8}

Ejemplo 1 Desarrollo con Sympy-Python

Con el algoritmo se tiene el siguiente resultado:

busca_PolosCeros:
 Q_polos : {-1: 1, -2 - 2*I: 1, -2 + 2*I: 1}
 P_ceros : {3: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 0

 asintoticamente:  estable

 h(t): 
     -t                   -2*t                            -2*t                
  4*e  *Heaviside(t)   9*e    *sin(2*t)*Heaviside(t)   4*e    *cos(2*t)*Heaviside(t)
- ------------------ + ----------------------------- + -----------------------------
          5                          10                              5        
>>>

En la gráfica se observa que todos los polos se encuentran en el lado izquierdo del plano LHP. Por lo que se considera que existen asintotas en la función correspondiente a tiempo, y el sistema es «asintóticamente estable«. También se considera BIBO estable por no tener los polos en el lado derecho RHP=0

Hs estabilidad polos Ej01

la gráfica de h(t) muestra lo indicado por el análisis de polos

Hs estabilidad ht Ej01

Instrucciones en Python Ejemplo 1

El algoritmo es una continuación del manejo de Polos y Ceros realizado en la sección anterior, donde se añade la revisión de las raíces de Q cuando son de tipo imaginarias.

Por simplicidad para el análisis se usan solo las partes esenciales, asi el enfoque es sobre la posición de los polos y su relación con h(t).

En el proceso, se agrupan las partes reales en Q_real y las imaginaria en Q_imag para realizar el conteo de lo repetido.

# Polos y ceros de funcion de transferencia H(s)
# Estabilidad del sistema, H(s) respuesta al impulso
# Ps es numerador, Qs es denominador
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm  

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
polosceros = fcnm.busca_polosceros(Hs)
P_ceros = polosceros['P_ceros']
Q_polos = polosceros['Q_polos']

# h(t) desde H_fp(s) con inversa Laplace 
ht = sym.inverse_laplace_transform(Hs,s,t)

# SALIDA

# Analiza estabilidad asintotica
casi_cero = 1e-12
# Separa parte real e imaginaria de raices
cuenta_real = 0; cuenta_imag = 0
unicos = 0 ; repetidos = 0 ; enRHP = 0
for raiz in Q_polos:
    [r_real,r_imag] = raiz.as_real_imag()
    if abs(r_real)>casi_cero and abs(r_imag)<casi_cero :
        cuenta_real = cuenta_real+1
    # para estabilidad asintotica
    conteo = Q_polos[raiz]
    if conteo==1 and r_real==0 and abs(r_imag)>0:
        unicos = unicos + 1
    if conteo>1  and r_real==0 and abs(r_imag)>0:
        repetidos = repetidos + 1
    if r_real>0:
        enRHP = enRHP + 1
cuenta_imag = len(Q_polos)-cuenta_real

# Revisa lado derecho del plano RHP
asintota = ""
if enRHP==0:
    asintota = 'estable'
if enRHP>0 or repetidos>0:
    asintota = 'inestable'
if enRHP==0 and unicos>0:
    asintota = 'marginalmente estable'

# SALIDA
print('busca_PolosCeros:')
for k in polosceros.keys():
    eq_lista = ['H','H_fp']
    if k in eq_lista:
        print('',k,':')
        sym.pprint(polosceros[k])
    else:
        print('\n',k,':',polosceros[k])

print(' polos reales: ',cuenta_real)
print(' polos complejos: ',cuenta_imag)
print(' sobre lado derecho RHP:',enRHP)
print(' sobre Eje Imaginario, repetidos: ',
      repetidos,' unicos:', unicos)
print('\n asintoticamente: ', asintota)

print('\n h(t): ')
sym.pprint(ht)

# GRAFICA -----------
muestras_H = 161
fig_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,
                          muestras=muestras_H,
                          f_nombre='H',
                          solopolos=True)
figura_h = fcnm.graficar_ft(ht,t_a,t_b,f_nombre='h')
plt.show()

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]
..


Ejemplo 2 H(s) con polos no repetidos, uno en RHP

Referencia: Lathi 2.14b p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D-1)(D^2 + 4D +8) y(t) = (D+2) x(t)

Ejemplo 2 Desarrollo analítico

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s-1)(s^2 + 4s +8) Y(s) = (s+2) X(s) \frac{Y(s)}{X(s)} = \frac{s+2}{(s-1)(s^2 + 4s +8)}

El grado del polinomio Q es N=3 y es mayor que el grado del polinomio P es M=1.

Los polos se obtienen como las raíces del denominador Q, necesario para plantear las fracciones parciales y simplificar la expresión en componentes más simples.

La expresión para el algoritmo se escribe como:

Hs = ((s+2)/(s-1))*(1/(s**2+4*s+8))

y el resultado de polos obtenido es:

polos: {1: 1, -2 - 2*I: 1, -2 + 2*I: 1}

Las fracciones parciales de H(s) usando el algoritmo de la sección anterior se obtiene

H(s) = \frac{3}{13}\frac{1}{s-1}+\frac{1}{13} \frac{3s+2}{s^2 + 4s +8}

Ejemplo 2 Desarrollo con Sympy-Python

Con el algoritmo se tiene el siguiente resultado:

busca_PolosCeros:
 H :
        s + 2         
----------------------
        / 2          \
(s - 1)*\s  + 4*s + 8/
 H_fp :
       3*s + 2            3     
- ----------------- + ----------
     / 2          \   13*(s - 1)
  13*\s  + 4*s + 8/             

 P_ceros : {-2: 1}

 Q_polos : {1: 1, -2 - 2*I: 1, -2 + 2*I: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 1
 sobre Eje Imaginario, repetidos:  0  unicos: 0

 asintoticamente:  inestable

 h(t): 
   t                   -2*t                            -2*t                   
3*e *Heaviside(t)   2*e    *sin(2*t)*Heaviside(t)   3*e    *cos(2*t)*Heaviside(t)
----------------- + ----------------------------- - -----------------------------
        13                        13                              13          
>>>

En la gráfica se observa que un polo se encuentra en el lado derecho del plano RHP. Un polo en RHP genera un componente creciente en el tiempo, en consecuencia, el sistema es «asintóticamente inestable». También el sistema es BIBO inestable.
Hs estabilidad polos Ej02

la gráfica de h(t) complementa lo interpretado con la posición de los polos en el plano s

Hs estabilidad ht Ej02

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]

..


Ejemplo 3 H(s) con polos complejos sobre eje imaginario únicos

Referencia: Lathi 2.14c p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D+2)(D^2 +4)y(t) = (D^2+D+1)x(t)

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+2)(s^2 + 4) Y(s) = (s^2+s+1) X(s)
Ps = s**2+s+1
Qs = (s+2)*(s**2+4)
Hs = Ps/Qs

El grado del polinomio Q es N=3 y el grado del polinomo P  es M=2, los polos del denominador muestran raíces complejas, por lo que se tendrá componentes cuadráticos para la transformada inversa de Laplace.

 Q_polos : {-2: 1, -2*I: 1, 2*I: 1}

Usando fracciones parciales con el algoritmo se tiene que,

H(s) = \frac{3}{8}\frac{1}{s+2}+\frac{1}{8} \frac{5s-2}{s^2 + 4}

resultados:

busca_PolosCeros:
 Q_polos : {-2: 1, -2*I: 1, 2*I: 1}
 P_ceros : {-1/2 - sqrt(3)*I/2: 1, -1/2 + sqrt(3)*I/2: 1}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 2

 asintoticamente:  marginalmente estable

 h(t): 
                                                       -2*t             
  sin(2*t)*Heaviside(t)   5*cos(2*t)*Heaviside(t)   3*e    *Heaviside(t)
- --------------------- + ----------------------- + --------------------
            8                        8                       8          
>>>

al tener polos sobre el eje imaginario el sistema es asintóticamente marginal estable, pero BIBO-inestable.

Hs estabilidad polos Ej03

con gráfica de h(t)

Hs estabilidad ht Ej03

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


Ejemplo 4. H(s) con polos repetidos sobre eje imaginario

Referencia: Lathi 2.14d p201

Revise la estabilidad BIBO de los siguientes sistemas LTIC descritos por las ecuaciones:

(D+1)(D^2 + 4)^2 y(t) = (D^2 + 2D + 8) x(t)

Se expresa la ecuación de operadores diferenciales D por la transformada de Laplace.

(s+1)(s^2 + 4)^2 Y(s) = (s^2 + 2s + 8) X(s)

Las raíces del denominador Q muestran que tienen valores complejos, mostrando que se usaría un componente cuadratico para la transformada inversa de Laplace.

 polos: {-1: 1, -2*I: 2, 2*I: 2}

al presentar dos raíces repetidas sobre el eje imaginario, se considera asintóticamente inestable y también BIBO inestable.

Revise el resultado de la inversa de Laplace usando las tablas.

los resultados con el algoritmo son:

busca_PolosCeros:
 H :
    2            
   s  + 2*s + 8  
-----------------
                2
        / 2    \ 
(s + 1)*\s  + 4/ 
 H_fp :
   2*(s - 6)     7*(s - 1)        7     
- ----------- - ----------- + ----------
            2      / 2    \   25*(s + 1)
    / 2    \    25*\s  + 4/             
  5*\s  + 4/                            

 P_ceros : {-1 - sqrt(7)*I: 1, -1 + sqrt(7)*I: 1}

 Q_polos : {-1: 1, -2*I: 2, 2*I: 2}
 polos reales:  1
 polos complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  2  unicos: 0

 asintoticamente:  inestable        
>>>

gráfica de polos

Hs estabilidad polos Ej04

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


Tarea

Lathi ejercicios 2.15 p202


Ejemplo 5. Con suma de términos y exponenciales de retraso en tiempo

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

la expresión para el algoritmo se escribe como:

Hs = (s+2)/(s-1) + sym.exp(-2*s)/(s-2)

con lo que se obtiene el resultado, usando intervalo de tiempo [-1,3] para las gráficas.

exp(-2*s) : {'Q_polos': {2: 1}, 'P_ceros': {}, 'Hs_k': 1/(s - 2)}
1 : {'Q_polos': {1: 1}, 'P_ceros': {-2: 1}, 'Hs_k': (s + 2)/(s - 1)}
Q_polos : {2: 1, 1: 1}
P_ceros : {-2: 1}
n_polos_real : 2
n_polos_imag : 0
enRHP : 2
unicos : 0
repetidos : 0
asintota : inestable

 h(t): 
 -4  2*t                       t                             
e  *e   *Heaviside(t - 2) + 3*e *Heaviside(t) + DiracDelta(t)
>>>

con gráfica de polos

Hs estabilidad polos Ej05

con gráfica h(t)

Hs estabilidad ht Ej05

Instrucciones con Python

# Estabilidad del sistema, H(s) respuesta al impulso
# Polos y ceros de funcion de transferencia H(s)
# Ps es numerador, Qs es denominador
import sympy as sym
import telg1001 as fcnm  

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

Hs = (s+2)/(s-1) + sym.exp(-2*s)/(s-2)

# Ps = s**2+2*s+8
# Qs = (s+1)*((s**2+4)**2)
# Hs = Ps/Qs

#Ps = s**2+s+1
#Qs = (s+2)*(s**2+4)
#Hs = Ps/Qs

#Hs = ((s+2)/(s-1))*(1/(s**2+4*s+8))
#Hs = ((s-3)/(s+1))*(1/(s**2+4*s+8))

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs_fp = fcnm.apart_s(Hs)
polosceros = fcnm.busca_polosceros(Hs_fp)
P_ceros = polosceros['P_ceros']
Q_polos = polosceros['Q_polos']

def estabilidad_asintotica_s(Q_polos, casi_cero=1e-8):
    ''' Analiza estabilidad asintotica con Q_raiz
        Separa parte real e imaginaria de raices
        casicero es la tolerancia para considerar cero
    '''
    # Separa parte real e imaginaria de raices
    cuenta_real = 0; cuenta_imag = 0
    unicos = 0 ; repetidos = 0 ; enRHP = 0
    for raiz in Q_polos:
        [r_real,r_imag] = raiz.as_real_imag()
        if abs(r_real)>casi_cero and abs(r_imag)<casi_cero :
            cuenta_real = cuenta_real+1
        # para estabilidad asintotica
        conteo = Q_polos[raiz]
        if conteo==1 and r_real==0 and abs(r_imag)>0:
            unicos = unicos + 1
        if conteo>1  and r_real==0 and abs(r_imag)>0:
            repetidos = repetidos + 1
        if r_real>0:
            enRHP = enRHP + 1
    cuenta_imag = len(Q_polos)-cuenta_real

    # Revisa lado derecho del plano RHP
    asintota = ""
    if enRHP==0:
        asintota = 'estable'
    if enRHP>0 or repetidos>0:
        asintota = 'inestable'
    if enRHP==0 and unicos>0:
        asintota = 'marginalmente estable'
        
    estable = {'n_polos_real': cuenta_real,
               'n_polos_imag': cuenta_imag,
               'enRHP'     : enRHP,
               'unicos'    : unicos,
               'repetidos' : repetidos,
               'asintota'  : asintota,}
    return(estable)

estable = estabilidad_asintotica_s(Q_polos)

# h(t) desde H_fp(s) con inversa Laplace
ht = 0*t
term_suma = sym.Add.make_args(Hs_fp.expand())
for term_k in term_suma:
    ht_k = sym.inverse_laplace_transform(term_k,s,t)
    ht = ht +ht_k
ht = sym.expand(ht,t) # terminos suma

# SALIDA
fcnm.print_resultado_dict(polosceros)
fcnm.print_resultado_dict(estable)

print('\n h(t): ')
sym.pprint(ht)

# GRAFICA -----------
import matplotlib.pyplot as plt
muestras_H = 161
fig_Hs = fcnm.graficar_Fs(Hs_fp,Q_polos,P_ceros,
                          muestras=muestras_H,
                          f_nombre='H',solopolos=True)
figura_h = fcnm.graficar_ft(ht,t_a,t_b,f_nombre='h')
plt.show()

polos: [ no repetidos con RHP=0]  [ no repetidos con RHP>=1 ]  [ complejos sobre eje Imag únicos ]  [ complejos sobre eje Imag repetidos ]  [ H(s) con exponencial ]


H(s) Implicaciones de estabilidad y Osciladores

Referencia: Lathi 2.5-3 p203

Todos los sistemas de procesamiento de señales deben ser asintóticamente estables.

Los sistemas inestables no son requeridos debido a que con condiciones iniciales, intencionales o no, llevan a respuestas fuera de rango. Estas respuestas inestables destruyen el sistema o lo llevan a condiciones de saturación que cambian la naturaleza del sistema debido al tipo de crecimiento exponencial.

Sistemas marginalmente estables, aunque BIBO inestables, tienen aplicación como osciladores, que es un sistema que genera su propia señal sin la aplicación de entradas externas.

La salida de un oscilador es una respuesta de entrada cero. Si la respuesta es una sinusiode con frecuencia ω0  el sistema es marginalmente estable en ±jω0. Para diseñar un oscilador de frecuencia  se debe usar un sistema descrito por la ecuación diferencial:

(D^2 + \omega _0 ^2 ) y(t) = x(t)

Sin embargo, osciladores en la práctica se construyen con sistemas no lineales.

4.3.2 H(s) – Polos reales y complejos con Sympy-Python

Referencia: Lathi ejemplo 4 p330, Hsu literal C. p112

La función de transferencia H(s) se escribe como P(s)/Q(s) que es una expresión racional de s. Se considera que en la expresión el grado m del polinomio del numerador P(s)  es mayor que el grado n del polinomio del denominador Q(s). Las raíces del polinomio P(s) del numerador se conocen como ceros, pues convierten la expresión en cero, mientras que para el denominador las raices se denominan polos al convertir en infinito la expresión (Hsu).

\frac{P(s)}{Q(s)} = \frac{a_0(s-ceros(1))\text{...}(s-ceros(m))}{b_0(s - polos(1)) \text{...} (s - polos(n))}

otra forma de escribir los polos y ceros (Lathi)

\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)

Para el desarrollo de un algoritmo con Sympy, el primer paso consiste en determinan los polos de un ejercicio.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo, busca_polosceros(Hs) ]


Ejemplo 1. Polos reales no repetidos y fracciones parciales

Referencia: Hsu problema 3.17a p137, Lathi ejemplo 4.3a p338, Oppenheim ejemplo 9.9 p671

La expresión H(s)  tiene los componentes de P(s) y Q(s) como se indican,

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

Ingreso de H(s) como P(s) y Q(s)

Para la expresión de la función de transferencia, con Sympy se crean los polinomios del numerador Ps y denominador Qs.

# INGRESO
s = sym.Symbol('s')
Ps = 2*s + 4  # 1+0*s cuando es constante
Qs = s**2 + 4*s + 3

Hs = Ps/Qs

Siendo H(s) la forma en que se expresan los ejercicios, se realizan las operaciones para obtener las raíces del numerador y denominador. El resultado buscado con el algoritmo es:

H(s):
   2*(s + 2)   
---------------
(s + 1)*(s + 3)

 H(s) en fracciones parciales 
  1       1  
----- + -----
s + 3   s + 1

 {Q_polos:veces}: {-1: 1, -3: 1}
 {P_ceros:veces}: {-2: 1}
>>>  

Instrucciones en Python

Hs se convierte en una expresión de Sympy con sym.sympify(Hs)s para el caso en que se de tan solo una constante. En el caso que Hs requiera agrupar términos se usa la instrucción sym.simplify(Hs,inverse=True).

Una vez preparada la expresión Hs, se separa como numerador Ps y denominador Qs en forma de polinomio para obtener las raíces y las veces que ocurren con la instrucción sym.roots().

La expresión de H(s) en fracciones parciales se obtienen usando sym.apart(Hs,s).

# Busca Polos y Ceros de H(s) con Sympy
# Ps es numerador, Qs es denominador
import sympy as sym

# INGRESO
s = sym.Symbol('s')

Ps = 2*s + 4
Qs = s**2 + 4*s + 3
Hs = Ps/Qs

# PROCEDIMIENTO
Hs = sym.sympify(Hs) # convierte a sympy una constante
Hs = sym.simplify(Hs,inverse=True) # agrupa terminos

# polos y ceros de Hs
[P,Q] = Hs.as_numer_denom()
P = P.as_poly(s)  # numerador
Q = Q.as_poly(s)  # denominador
P_ceros = sym.roots(P)
Q_polos = sym.roots(Q)

# en factores
Hs = sym.factor(Hs,s)
# fracciones parciales
Hs_fp = sym.apart(Hs,s)

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en fracciones parciales ')
sym.pprint(Hs_fp)
print('\n {Q_polos:veces}:',Q_polos)
print('\n {P_ceros:veces}:',P_ceros)

Para visualizar el concepto de polos y ceros, se presenta la gráfica de polos y ceros en el dominio s, para un corte en el plano real. Por simplicidad se usará el procedimiento creado para el curso en telg1001.py

Todos los polos y ceros se encuentran en el lado izquierdo del plano, es decir su parte real es negativa.

# GRAFICA -----------
import matplotlib.pyplot as plt
import telg1001 as fcnm
fig_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,f_nombre='H')
plt.show()

Hs Polos Ceros Ej01

Para revisar los polos de la expresión de fracciones parciales, se considera que los términos de las expresiones podrían tener polos repetidos como en el ejemplo 3. Para polos repetidos, se usa el número de veces mayor cuando aparecen polos de términos suma.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]

..


Ejemplo 2. Grados iguales de P(s) y Q(s), ceros con valores complejos

Referencia: Lathi ejemplo 4.3b p338, Hsu ejercicio 3.20 p140, Oppenheim ejemplo 9.36 p716

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

Ejemplo 2 Desarrollo con Sympy-Python

Al algoritmo del ejercicio anterior se modifica solo el bloque de ingreso con las expresiones para numerador y denominador:

Ps = 2*s**2+5
Qs = s**2 + 3*s + 2
Hs = Ps/Qs

El resultado Hs_fp en fracciones parciales muestra que existe un término constante en las sumas de la expresión. La forma descrita por Lathi al inicio de la página, la constante se la considera como ‘ganancia(s)’.

H(s):
       2       
    2*s  + 5   
---------------
(s + 1)*(s + 2)

 H(s) en fracciones parciales 
      13      7  
2 - ----- + -----
    s + 2   s + 1


 {Q_polos:veces}: {-1: 1, -2: 1}
 {P_ceros:veces}: {-sqrt(10)*I/2: 1, sqrt(10)*I/2: 1}
>>>  

Por otra parte se tiene que los ceros tienen raíces con componente imaginario representado con el símbolo ‘I‘, que para representar cada raíz en la gráfica se considera el eje vertical también como parte Imag(s). Si es cierto, es diferente a la naturaleza de H(s), perdonen, se usa ésta «libertad» solo para resaltar la influencia de los polos en H(s) como se mostró en el ejercicio anterior, o se realizaría un gráfico con el plano real e imaginario como base y una tercera dimensión para H(s).

Los polos se encuentran en el lado izquierdo del plano, es decir su parte real es negativa.

Hs Polos Ceros Ej02

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]


Ejemplo 3. Polos reales repetidos o raices repetidas en denominador Q(s)

Referencia: Lathi ejemplo 4.3d p342

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

la forma de la expresión para el algoritmo es

Ps = 8*s+10
Qs = (s+1)*((s+2)**3)

Hs = Ps/Qs

En él ejercicio el término del denominador tiene un factor elevado al cubo, lo que muestra que la raiz es repetida y se debe reflejar en la respuesta del algoritmo. En las fracciones parciales se muestra que para el polo repetido existe un término suma para cada grado en el denominador.

H(s):
  2*(4*s + 5)   
----------------
               3
(s + 1)*(s + 2) 

 H(s) en fracciones parciales 
    2        2          6         2  
- ----- - -------- + -------- + -----
  s + 2          2          3   s + 1
          (s + 2)    (s + 2)         

 {Q_polos:veces}: {-1: 1, -2: 3}
 {P_ceros:veces}: {-5/4: 1}
>>> 

Todos los polos y ceros se encuentran en el lado izquierdo del plano s.

Se adjunta la gráfica para la interpretación de H(s)

H(s) Polos Ceros Ej03

Si la busqueda de polos es sobre la expresión de fracciones parciales, se tiene que las veces se determinan como el mayor exponente del término del denominador.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]

..


Ejemplo 4. H(s) con Fracciones parciales y polos de tipo complejo

Referencia: Lathi ejemplo 4.3c p340

H(s) = \frac{6(s+34)}{s(s^2+10s+34)}

la forma de la expresión para el algoritmo es

Ps = 6*(s+34)
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

Incorporando las instrucciones para identificar si el denominador tiene raices complejas, se asignan los valores de los parametros necesarios en un diccionario de variables, obteniendo como resultado,

H(s):
    6*(s + 34)    
------------------
  / 2            \
s*\s  + 10*s + 34/

 H(s) en fracciones parciales 
    6*(s + 9)      6
- -------------- + -
   2               s
  s  + 10*s + 34    

 {Q_polos:veces}: {-5 - 3*I: 1, -5 + 3*I: 1, 0: 1}
 {P_ceros:veces}: {-34: 1}

 h(t): 
                               -5*t                              
- 2*(4*sin(3*t) + 3*cos(3*t))*e    *Heaviside(t) + 6*Heaviside(t)
>>> 

grafica ejemplo 4

H(s) Polos Ceros Ej04

Ejemplo 4 Instrucciones en Python

# Fracciones parciales de H(s) con Sympy
# Ps es numerador, Qs es denominador
# con Transformada Inversa de Laplace
import numpy as np
import sympy as sym

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

Ps = 6*(s+34) 
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

# PROCEDIMIENTO
Hs = sym.simplify(Hs,inverse=True)
# Hs debería ser un solo termino suma
Hs = sym.factor(Hs,s)
# fracciones parciales
Hsp = sym.apart(Hs,s)

# Analiza polos
def polosceros_simple(Hs):
    ''' Busca_polo de un termino sin exp(-a*s)
    '''
    Hs = sym.simplify(Hs,inverse=True)
    # Hs debería ser un solo termino suma
    Hs_factor = sym.factor(Hs,s)
    # polos y ceros de termino Hs
    [P,Q] = Hs_factor.as_numer_denom()
    P = sym.poly(P,s)  # numerador
    Q = sym.poly(Q,s)  # denominador
    P_ceros = sym.roots(P)
    Q_polos = sym.roots(Q)
    respuesta = {'Q_polos' : Q_polos,
                 'P_ceros' : P_ceros}
    return(respuesta)

polosceros = polosceros_simple(Hs)

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en fracciones parciales ')
sym.pprint(Hsp)
print('\n {Q_polos:veces}:',polosceros['Q_polos'])
print(' {P_ceros:veces}:',polosceros['P_ceros'])
print('\n h(t): ')
sym.pprint(ht)

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]


Ejemplo 5. Con suma de términos y exponenciales de retraso en tiempo

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

la expresión para el algoritmo se escribe como:

Hs = ((s+2)/(s-1)) + sym.exp(-2*s)/(s-2)

Para el segundo término que tiene un exponencial, no se puede aplicar directamente la instrucción sym.apart(Hs,s) al no ser reconocida como tipo polinomio. Para continuar se debe separar el término sym.exp() del resto de la expresión, y se podrá aplicar fracciones parciales a esa parte, para luego volver a poner el término.

Para facilitar el desarrollo del algoritmo, se recomienda empezar a analizar una expresión mas simple, como solamente el segundo término de la suma.

 expresion H(s):
         -2*s
s + 2   e    
----- + -----
s - 1   s - 2

Polos y ceros:
término: 1 * (s + 2)/(s - 1)
 Q_polos{polos:veces}:  {1: 1}
 P_ceros{polos:veces}:  {-2: 1}
término: exp(-2*s) * 1/(s - 2)
 Q_polos{polos:veces}:  {2: 1}
 P_ceros{polos:veces}:  {}
>>> 

H(s) Polos Ceros Ej05

En el algoritmo se empieza por determinar si la expresión de Hs tiene términos sym.exp(-s), la instrucción que permite listar estos términos es

lista_exp = list(Hs.atoms(sym.exp(-s)))

como la expresión H(s) en los ejercicios puede contener sumas o multiplicaciones de bloques, primero se simplifica la expresión H(s)

    # convierte a Sympy si es solo constante
    Hs = sym.sympify(Hs)
    # simplifica operación Hs1(s)*Hs2(s)
    Hs = sym.expand_power_exp(Hs)
    Hs = sym.expand(Hs, power_exp=False)
    
    term_sum = sym.Add.make_args(Hs)
    Hs_0 = 0*s # reagrupa Fs
    for term_k in term_sum:
        term_k = sym.simplify(term_k)
        term_k = sym.expand_power_exp(term_k)
        Hs_0 = Hs_0 + term_k
    
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_0.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

luego se tabula la lista_exp y se eliminan los términos que pueden ser constantes o no dependen de ‘s’, por ejemplo sym.exp(10).

Luego se separa la expresión por grupos usando sym.collect() creando un diccionario con los elementos de la lista_exp y con las expresiones agrupadas.

    # separados por lista_exp
    separados = sym.collect(Hs_0,lista_exp,evaluate=False)

Con los elementos separados se puede aplicar la búsqueda de polos y ceros para cada elemento de la lista con la función polosceros_simple(Hs_k)

Al final se pueden agrupar los Q_polos y P_ceros de toda la expresión siguiendo criterios como que las veces que ocurre cada polo es el maximo entre los términos separados por lista_exp.

Instrucciones en Python

# Busca Polos y Ceros de H(s) con Sympy
# separa términos con sym.exp(-a*s)
# Ps es numerador, Qs es denominador
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')

# Ps = (2*s + 4)
# Qs = s**2 + 4*s + 3

# Ps = 2*s**2+5
# Qs = s**2 + 3*s + 2

# Ps = 8*s+10
# Qs = (s+1)*((s+2)**3)

# Ps = 6*(s+34)
# Qs = s*(s**2+10*s+34)

# Hs = Ps/Qs
# Hs = sym.exp(-2*s)/(s-2)
# Hs = 1/s - sym.exp(-2*s)/s
# Hs = sym.exp(-2*s)/s -sym.exp(-4*s)/s
# Hs = (2*s**2+5*s+12)/((s+2)*(s**2+2*s+10))
Hs = ((s+2)/(s-1)) + sym.exp(-2*s)/(s-2)

# PROCEDIMIENTO
def busca_polosceros(Hs):
    ''' Busca polos de Hs (divisiones para cero)
        y ceros de Hs,cuenta las veces que aparecen,
        agrupa por exp(-a*s) agrupado en factores.
    '''
    def polosceros_simple(Hs):
        ''' Busca_polo de un termino sin exp(-a*s)
        '''
        Hs = sym.simplify(Hs,inverse=True)
        # Hs debería ser un solo termino suma
        Hs_factor = sym.factor(Hs,s)
        # polos y ceros de termino Hs
        [P,Q] = Hs_factor.as_numer_denom()
        P = sym.poly(P,s)  # numerador
        Q = sym.poly(Q,s)  # denominador
        P_ceros = sym.roots(P)
        Q_polos = sym.roots(Q)
        respuesta = {'Q_polos' : Q_polos,
                     'P_ceros' : P_ceros}
        return(respuesta)
    
    # convierte a Sympy si es solo constante
    Hs = sym.sympify(Hs)
    # simplifica operación Hs1(s)*Hs2(s)
    Hs = sym.expand_power_exp(Hs)
    Hs = sym.expand(Hs, power_exp=False)
    
    term_sum = sym.Add.make_args(Hs)
    Hs_0 = 0*s # reagrupa Fs
    for term_k in term_sum:
        term_k = sym.simplify(term_k)
        term_k = sym.expand_power_exp(term_k)
        Hs_0 = Hs_0 + term_k
    
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_0.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

    # separados por lista_exp
    separados = sym.collect(Hs_0,lista_exp,evaluate=False)

    # polos y ceros por terminos exp(-s) agrupados
    polosceros = {} ; Hs_fp = 0
    for k in separados:
        Hs_k = sym.factor(separados[k],s)
        PZ_k = polosceros_simple(Hs_k)
        PZ_k['Hs_k']  = Hs_k
        polosceros[k] = PZ_k

    # integra polos
    Q_polos = {} ; P_ceros = {}
    for k in polosceros:
        polos = polosceros[k]['Q_polos']
        for unpolo in polos:
            if unpolo in Q_polos:
                veces = max([Q_polos[unpolo],polos[unpolo]])
                Q_polos[unpolo] = veces
            else:
                Q_polos[unpolo] = polos[unpolo]
        ceros = polosceros[k]['P_ceros']
        for uncero in ceros:
            if uncero in P_ceros:
                veces = max([P_ceros[uncero],ceros[uncero]])
                P_ceros[unpolo] = veces
            else:
                P_ceros[uncero] = ceros[uncero]
    # revisa exp(-a*s)
    if len(lista_exp)==0: #sin componentes
        del polosceros['1']
    polosceros['Q_polos'] = Q_polos
    polosceros['P_ceros'] = P_ceros
    return(polosceros)

polosceros = busca_polosceros(Hs)

# SALIDA
print('\n expresion H(s):')
sym.pprint(Hs)
print('\nPolos y ceros:')
lista_PZ = ['Q_polos','P_ceros']
for k in polosceros:
    if not(k in lista_PZ):
        print('término:', k,'*',polosceros[k]['Hs_k'])
        print(' Q_polos{polos:veces}: ',polosceros[k]['Q_polos'])
        print(' P_ceros{polos:veces}: ',polosceros[k]['P_ceros'])
    else:
        print(k,polosceros[k])

La gráfica en éste caso se realiza usando Hs y sus componentes, de ésta forma es posible visualizar el efeto de cada componente con sus polos y ceros.

# GRAFICA ----------------------
def graficar_Fs2(Fs,Q_polos={},P_ceros={},
                s_a=1,s_b=0,muestras=101,f_nombre='F',
                solopolos=False,polosceros={}):
    # grafica Fs
    lista_Hs =[]

    lista_PZ = ['Q_polos','P_ceros']
    lista_term = list(polosceros.keys())
    lista_term.remove('Q_polos')
    lista_term.remove('P_ceros')

    # intervalo para s,y componentes Fs con exp(-a*s)
    s_a = 0 ; s_b = 0 ; Lista_Hs = []
    cond_componente = False
    if len(lista_term)>0:
        cond_componente = True
        if len(polosceros)>1: # tiene componenes con exp(-a*s)
            lista_Hs.append([Hs,{},{}])
        for k in polosceros:
            if not(k in lista_PZ):
                Q_polosk = polosceros[k]['Q_polos']
                P_cerosk = polosceros[k]['P_ceros']
                [s_ak,s_bk] = fcnm.intervalo_s(Q_polosk,P_cerosk)
                s_a = min(s_a,s_ak)
                s_b = max(s_b,s_bk)
                lista_Hs.append([polosceros[k]['Hs_k']*sym.sympify(k),
                                 Q_polosk,P_cerosk])
    else: #sin componenes con exp(-a*s)
        lista_Hs.append([Hs,polosceros['Q_polos'],polosceros['P_ceros']])
        [s_a,s_b] = fcnm.intervalo_s(polosceros['Q_polos'],
                                       polosceros['P_ceros'])
    # graficas por cada componente
    fig_Fs, graf_Fs = plt.subplots()
    # no presenta errores de división para cero en lambdify()
    np.seterr(divide='ignore', invalid='ignore')
    for k in lista_Hs:
        Fs_k = k[0] ;Q_polosk = k[1] ;P_cerosk = k[2]
        
        # convierte a sympy una constante
        Fs_k = sym.sympify(Fs_k,s) 
        if Fs_k.has(s): # no es constante
            F_s = sym.lambdify(s,Fs_k,modules=fcnm.equivalentes)
        else:
            F_s = lambda s: Fs_k + 0*s
        
        s_i = np.linspace(s_a,s_b,muestras)
        Fsi = F_s(s_i) # Revisar cuando s es complejo
        lineaestilo = 'solid'
        if len(Q_polosk)>0 and cond_componente:
            lineaestilo = 'dashed'
            
        graf_Fs.plot(s_i,Fsi,label=Fs_k,linestyle=lineaestilo)
        lineacolor = plt.gca().lines[-1].get_color()

        if len(Q_polosk)>0:
            polo_re = [] ; polo_im = []
            for polos in Q_polosk.keys():    
                graf_Fs.axvline(sym.re(polos),color='red',
                                linestyle='dotted')
                x_polo = sym.re(polos)
                y_polo = sym.im(polos)
                polo_re.append(x_polo)
                polo_im.append(y_polo)
                etiqueta = "("+str(x_polo)+','+str(y_polo)+")"
                graf_Fs.annotate(etiqueta,(x_polo,y_polo))
        
            graf_Fs.scatter(polo_re,polo_im,marker='x',
                            color =lineacolor,
                            label = Q_polosk)
        if len(P_cerosk)>0:
            cero_re = [] ; cero_im = []
            for cero in P_cerosk.keys():    
                x_cero = sym.re(cero)
                y_cero = sym.im(cero)
                cero_re.append(x_cero)
                cero_im.append(y_cero)
                etiqueta = "("+str(x_cero)+','+str(y_cero)+")"
                graf_Fs.annotate(etiqueta,(x_cero,y_cero))
        
            graf_Fs.scatter(cero_re,cero_im,marker='o',
                            color =lineacolor,
                            label = P_cerosk)

    graf_Fs.axvline(0,color='gray')
    graf_Fs.legend()
    graf_Fs.set_xlabel('Real(s)')
    graf_Fs.set_ylabel('F(s) ; Imag(s)')
    graf_Fs.set_title(r'F(s) = $'+str(sym.latex(Hs))+'$')
    graf_Fs.grid()
    return(fig_Fs)

fig_Fs = graficar_Fs2(Hs,polosceros=polosceros)
plt.show()

Parametros cuadráticos

Para el ejercicio 4, se tienen polos con valores complejos y para usar la tabla de transformadas de Laplace se requiere obtener algunos parámetros.

Referencia: Lathi ejemplo 4.3c p340

H(s) = \frac{6(s+34)}{s(s^2+10s+34)}

el desarrollo se puede describir como

H(s) = \frac{6(s+34)}{s(s^2+10s+34)} = \frac{k_1}{s} + \frac{As+B}{s^2+10s+34}

usando el método de Heaviside, se obtiene k1

k_1 = \frac{6(s+34)}{\cancel{s}(s^2+10s+34)}\Big|_{s=0} =\frac{6(0+34)}{(0^2+10(0)+34)} = \frac{6(34)}{34}=6 \frac{6(s+34)}{s(s^2+10s+34)} = \frac{6}{s} + \frac{As+B}{s^2+10s+34}

se simplifican las fracciones al multiplicar ambos lados por s(s2+10s+34)

6(s+34) = 6(s^2+10s+34) + (As+B)s 6s+204 = 6s^2+60s+204 + A s^2+Bs 6s+204 = (6+A)s^2+(60+B)s+204

se puede observar que 6+A=0, por lo que A=-6.
Tambien se tiene que 60+B=6, por lo que B=-54
quedando la expresión como:

H(s) = \frac{6}{s} + \frac{-6s-54}{s^2+10s+34}

Observando la tabla de transformadas de Laplace se compara la expresión con las filas 12c y 12d se usan los parametros de:

\frac{As+B}{s^2+2as+c}

A= -6, B=-64, a=5, c = 34
y se calculan:

b = \sqrt{c-a^2} = \sqrt{34-(5)^2} = 3 r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} r = \sqrt{\frac{(-6)^2 (34) +(-64)^2 -2(-6)(-64)(5)}{(34)-(5)^2}}=10 \theta = \tan ^{-1} \Bigg( \frac{(-6)(5)-(-64)}{(-6)\sqrt{34-(5)^2}} \Bigg) = -0.9272

con lo que se puede escribir la transformada inversa de Laplace para h(t) como

h(t) = [6+10e^{-5t}\cos (3t--0.9272)] \mu (t)

o la otra foma de expresión usando la fila 12d de la tabla

Ejemplo 4 Desarrollo con Sympy-Python

Incorporando las instrucciones para identificar si el denominador tiene raices complejas, se asignan los valores de los parámetros necesarios en un diccionario de variables, obteniendo como resultado,

H(s) en fracciones parciales:
    6*(s + 9)      6
- -------------- + -
   2               s
  s  + 10*s + 34    

 parametros de Q_cuadratico: 
  -6*(s + 9)/(s**2 + 10*s + 34)  :
   {'A': -6.0, 'B': -54.0, 'a': 5.0, 'c': 34.0, 'r': 10.0, 'b': 3.0, 'theta': -0.9272952180016122}
>>> 

mientras que el algoritmo se resume en

# Fracciones parciales de H(s) con Sympy
# Ps es numerador, Qs es denominador
# con Transformada Inversa de Laplace
import numpy as np
import sympy as sym
import telg1001 as fcnm

# INGRESO
s = sym.Symbol('s')
t = sym.Symbol('t', real=True)

Ps = 6*(s+34) 
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

# PROCEDIMIENTO
# parametros de termino con Q cuadratico
def Q_cuad_s_parametros(Hs):
    '''Si Q tiene grado=2, obtiene los parametros
       para la tabla de Transformadas Inversas de Laplace
    '''
    def Q_cuad_sin_exp(Hs):
        '''para Hs sin sym.exp(-a*s)
        '''
        # fracciones parciales NO separa Q con factores complejos
        Hs_fp = sym.apart(Hs,s)
        respuesta = {}
        term_suma = sym.Add.make_args(Hs_fp)
        for term_k in term_suma:
            [P,Q]  = term_k.as_numer_denom()
            if sym.degree(Q) == 2 and sym.degree(P) == 1:
                P_coef = P.as_coefficients_dict(s)
                Q_coef = Q.as_coefficients_dict(s)
                Q0 = 1 # Normalizar coefficiente s**2=1
                if Q_coef[s**2]!=1: 
                    Q0 = Q_coef[s**2]
                # Parametros de Q cuadratico
                a = 0 ; c = 0
                if s**1 in Q_coef:
                    a = float(Q_coef[s**1]/Q0)/2
                if s**0 in Q_coef:
                    c = float(Q_coef[s**0]/Q0)
                A = float(P_coef[s**1])
                B = 0
                if s**0 in P_coef:
                    B = float(P_coef[s**0])
                rP = (A**2)*c + B**2 - 2*A*B*a
                rQ = c - a**2
                r  = np.sqrt(rP/rQ)
                b  = np.sqrt(c-a**2)
                thetaP = A*a-B
                thetaQ = A*np.sqrt(c-a**2)
                theta  = np.arctan(thetaP/thetaQ)
                parametro = {'A': A, 'B': B,'a': a, 'c': c,
                             'r': r, 'b': b,'theta':theta}
                respuesta[term_k] = parametro
        return (respuesta)

    Hs = sym.sympify(Hs,s)
    Hs_fp = fcnm.apart_s(Hs) # separa por sym.exp(-a*s)
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_fp.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

    # separados por lista_exp
    separados = sym.collect(Hs_fp,lista_exp,evaluate=False)

    Qs2 = {} # agrupa parametros Qs2
    for k in separados:
        Qs2_k = Q_cuad_sin_exp(separados[k])
        if len(Qs2_k)>0:
            for term in Qs2_k:
                Qs2[k*term] = Qs2_k[term]
    return(Qs2)

Hs_fp = fcnm.apart_s(Hs)
Qs2 = Q_cuad_s_parametros(Hs)

# SALIDA
print('H(s) en fracciones parciales:')
sym.pprint(Hs_fp)

# parametros de termino con Q cuadratico
if len(Qs2)>0:
    print('\n parametros de Q_cuadratico: ')
    for k in Qs2:
        print(' ',k,' :')
        print('  ',Qs2[k])

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]