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.
la gráfica de H(s) junto a los polos:
la gráfica de h(t) como respuesta al impulso,
la gráfica de la señal de entrada x(t) y la respuesta y(t) es:
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 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()