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

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

Se revisan los polos y ceros para escribir la función de transferencia,

H(s) = \frac{(s-ceros_0)(s-ceros_1)}{(s-polo_0)(s-polo_1)} H(s) = \frac{(s-(-1.5j))(s-(1.5j))}{(s-(-1+0.5j))(s-(-1-0.5j))}

H(s) = \frac{s^2+2.25}{s^2+2s+1.25}

tomando el modelo proporcionado de la ecuación, se compara para encontrar los valores de las variables

H(s) = \frac{k (s^2 + b_1 s + b_2)}{s^2 + a_1 s + a_2}

k = 1, b1=0, b2=2.25, a1=2, a2=1.25

Separando en fracciones parciales, dado que el grado M del numerador y grado N denominador son iguales, ganancia el coeficiente del término de mayor grado del numerador. Lo mismo que se obtiene al dividir el numerador y denominador por s2 y hacer s→∞

H(s) = \frac{1+\frac{2.25}{s^2}}{1+2\frac{1}{s}+\frac{1.25}{s^2}} \Big|_{s=\infty} =1 H(s) = 1+\frac{k_1s+k_2}{s^2+2s+1.25} H(s) = \frac{s^2+2s+1.25+k_1s+k_2}{s^2+2s+1.25} H(s) = \frac{s^2+(2+k_1)s+(1.25+k_2)}{s^2+2s+1.25}

con lo que 2+k1=0, entonces k1=-2,
también 1.25+k2=2.25, entonces K2 = 1

con lo que se obtiene,

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

los polos se encuentran sobre el lado izquierdo del plano s, por lo que el sistema es asintoticamente estable

Observe que para H(s) el grado del polinomio el numerador P de H(s) es igual al grado del denominador Q. por lo que M=N. por lo que se considera un función impropia (Lathi 4.3b p339, 4.2 p359). Se genera una respuesta impulso por el valor de la constante 1 que se obtiene al usar fracciones parciales.

Para la transformada inversa, para el término cuadrático, en la tabla de transformadas de Laplace se identifica el caso 12c:

\frac{As+B}{s^2+2as+c} donde A=-2, B=1,a=1,c=1.25

El valor b = \sqrt{c-a^2} se calcula como:

b = \sqrt{1.25-(1)^2} = \frac{1}{2} r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} r = \sqrt{\frac{(-2)^2 (1.25) +(1)^2 -2(-2)(1)(1)}{1.25-(-2)^2}} =\sqrt{\frac{10}{1/4}} = 6.3245 \theta = \tan ^{-1} \Big( \frac{Aa-B}{A\sqrt{c-a^2}}\Big) \theta = \tan ^{-1} \Big( \frac{(-2)(1)-(1)}{(-2)\sqrt{1.25-(1)^2}}\Big)\tan ^{-1}(\frac{-3}{-2(1/2)} = 1.2490

h(t) = δ(t) + re-atcos (bt+θ) μ(t)

h(t) = δ(t) + 6.3245e-tcos (t/2+1.2490) μ(t)

Algoritmo en Python

Para el bloque de ingreso se debe mantener los coeficientes como números enteros o racionales. Recordando que los algoritmos de Sympy para transformadas de Laplace (hasta versión 1.11) se encuentran escritos para el dominio de los enteros y racionales (‘ZZ’), por lo que los datos se ingresan como:

k1 = sym.Rational(2.25).limit_denominator(100)
k2 = sym.Rational(1.25).limit_denominator(100)
Hs = (s**2+k1)/(s**2+2*s+k2)

con el resultado mostrado,

 H(s) = P(s)/Q(s)
    2   9   
   s  + -   
        4   
------------
 2         5
s  + 2*s + -
           4

 H(s) fracciones parciales:
   4*(2*s - 1)      
- -------------- + 1
     2              
  4*s  + 8*s + 5    

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

 {polos,veces}:  {-1 - I/2: 1, -1 + I/2: 1}
 polos reales:  0  complejos:  2
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 0
 asintoticamente:  estable

 parametros Q cuadraticos s: 
{'A': -2.0, 'B': 1.0, 'a': 1.0, 'c': 1.25,
 'r': 6.324555320336759, 'b': 0.5,
 'theta': 1.2490457723982544}

Con gráficas de H(s) y polos

s1Eva2012TI_T1 polos H(s)

gráfica de h(t)

s1Eva2012TI_T1 h(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) y estabilidad
k1 = sym.Rational(2.25).limit_denominator(100)
k2 = sym.Rational(1.25).limit_denominator(100)
Hs = (s**2+k1)/(s**2+2*s+k2)
#Hs = (s**2+2.25)/(s**2+2*s+1.25)

# 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 = 5
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(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)

# 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])

# Graficas polos, H(s), con polos h(t) --------
muestras_H = 101
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,f_nombre='H')
figura_ht = fcnm.graficar_ft(ht,t_a,t_b,muestras,f_nombre='h')
plt.show()

.

s1Eva2009TII_T5 LTI DT bloques H[z] en serie

Referencia: 1Eva2009TII_T5 LTI DT bloques H[z] en serie

1. Las respuestas impulso de cada subsistema

H_1[z] = \frac{z}{z-(0.7)} = \frac{z}{z-0.7}

usando la tabla de transformadas z

h_1[n] = (0.7)^n \mu[n]

Continuando con el subsistema de la derecha

H_2[z] = \frac{z}{z-(-0.5)} = \frac{z}{z+0.5} h_2[n] = (-0.5)^n \mu[n]

El sistema total:

H[z] = H_1[z] H_2[z] = \frac{z}{(z-0.7)} \frac{z}{(z+0.5)} = \frac{z^2}{(z-0.7)(z+0.5)}

fracciones parciales modificadas, multiplica ambos lados por 1/z

\frac{H[z]}{z} = \Big( \frac{1}{z} \Big) \frac{z^2}{(z-0.7)(z+0.5)}= \frac{z}{(z-0.7)(z+0.5)} \frac{H[z]}{z} = \frac{z}{(z-0.7)(z+0.5)} = \frac{k_1}{z-0.7} +\frac{k_2}{z+0.5} k_1 = \frac{z}{\cancel{(z-0.7)}(z+0.5)} \Big|_{z=0.7} = \frac{0.7}{(0.7+0.5)} = 0.5833 k_2 = \frac{z}{(z-0.7)\cancel{(z+0.5)}} \Big|_{z=-0.5} = \frac{-0.5}{(-0.5-0.7)} = 0.4166 \frac{H[z]}{z} = \frac{0.5833}{z-0.7} +\frac{0.4166}{z+0.5}

Restaura fracciones parciales, multiplica ambos lados por z

H[z] = \frac{0.5833 z}{z-0.7} +\frac{0.4166z}{z+0.5}

usando la tabla de transformadas z

h[n] = 0.5833 \Big(0.7 \Big)^n \mu[n] +0.4166 \Big(-0.5\Big)^n \mu[n]

factor común μ[n]

h[n] = \Bigg( 0.5833 \Big(0.7 \Big)^n +0.4166 \Big(-0.5\Big)^n \Bigg) \mu[n]

Revisando el resultado con el algoritmo en Python

 Hz:
          2        
         z         
-------------------
(z - 0.7)*(z + 0.5)

 Hz/z:
          z          
---------------------
     2               
1.0*z  - 0.2*z - 0.35

 Hz/z.apart:
0.416666666666667   0.583333333333333
----------------- + -----------------
   1.0*z + 0.5         1.0*z - 0.7   

 Hz = (Hz/z)*z
0.416666666666667*z   0.583333333333333*z
------------------- + -------------------
    1.0*z + 0.5           1.0*z - 0.7    
 
polos:     {-0.500000000000000: 1, 0.700000000000000: 1}

 polos Re:  [-1/2, 7/10]

Tarea: 2. Su respuesta y[n]=s]n], expresada a la mínima expresión frente a la siguiente excitación x[n]=μ[n], esquematícela.

Instrucciones en Python

Algoritmos desarrollados en X[z] Fracciones parciales modificadas con Python

y la parte gráfica de Transformada z con Sympy-Python

# Transformada z- Fracciones parciales
# Polos únicos, repetidos y complejos
# Lathi Ejercicio 5.3c pd495
# blog.espol.edu.ec/telg1001
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

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

Pz = z**2
Qz = (z-0.7)*(z+0.5)

# PROCEDIMIENTO
P = Pz.as_poly(z)
Q = Qz.as_poly(z)
Xz = Pz/Qz

# Raices Denominador
Q_raiz = sym.roots(Q)
Q_raizRe = sym.real_roots(Q)
rcompleja = len(Q_raiz)-len(Q_raizRe)

# Raices reales, únicas y repetidas
if (rcompleja<=0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z

def parametro_cuadratico(untermino):
    unparametro ={}
    # revisa denominador cuadratico
    [numerador,denominador] = (untermino).as_numer_denom()
    gradoD = 0
    coeficientesD = denominador
    gradoN = 0
    coeficientesN = numerador
    if not(denominador.is_constant()):
        denominador = denominador.as_poly()
        gradoD = denominador.degree()
        coeficientesD = denominador.coeffs()
    if not(numerador.is_constant()):
        numerador = numerador.as_poly()
        gradoN = numerador.degree()
        coeficientesN = numerador.coeffs()
    if gradoD == 2 and gradoN==2:
        a = float(coeficientesD[1])/2
        gamma2 = float(coeficientesD[2])
        gamma = np.sqrt(gamma2)
        A = float(coeficientesN[0])
        B = float(coeficientesN[1])
        rN = (A**2)*gamma2 + B**2 - 2*A*a*B
        rD = gamma2 - a**2
        r = np.sqrt(rN/rD)
        beta = np.arccos(-a/gamma)
        thetaN = A*a-B
        thetaD = A*np.sqrt(gamma2-a**2)
        theta = np.arctan(thetaN/thetaD)
        unparametro = {'r':r,
                       'gamma':gamma,
                       'beta':beta,
                       'theta':theta}
    return (unparametro)

# Terminos cuadraticos
parametros = [] # denominador cuadratico
if (rcompleja>0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z
        unparam = parametro_cuadratico(untermino*z)
        if len(unparam)>0:
            parametros.append(unparam)
if (rcompleja>0 and len(Q_raizRe)==0):
    Xzp = P/Q
    parametros.append(parametro_cuadratico(P/Q))

# SALIDA
print('\n Hz:')
sym.pprint(Xz)
if (len(Q_raizRe)>0):
    print('\n Hz/z:')
    sym.pprint(Xzz)
    print('\n Hz/z.apart:')
    sym.pprint(Xzm)
    print('\n Hz = (Hz/z)*z')
    sym.pprint(Xzp)
    print('\n polos:    ', Q_raiz)
print('\n polos Re: ', Q_raizRe)
if len(parametros)>0:
    for unparam in parametros:
        print('parametros cuadraticos: ')
        for cadauno in unparam.keys():
            print(cadauno,'\t',unparam[cadauno])


# grafica plano z imaginario
figura, grafROC = plt.subplots()
# limite
radio1 = plt.Circle((0,0),1,color='lightgreen',
                    fill=True)
radio2 = plt.Circle((0,0),1,linestyle='dashed',
                    color='lightgreen',fill=False)
grafROC.add_patch(radio1)
for raiz in Q_raiz.keys():
    [r_real,r_imag] = raiz.as_real_imag()
    radio_raiz = np.abs(raiz)
    nROC = plt.Circle((0,0),radio_raiz,
                      color='red',fill=True)
    grafROC.add_patch(nROC)
grafROC.add_patch(radio2) # borde r=1
grafROC.axis('equal')
# marcas de r=1 y valor a
for raiz in Q_raiz.keys():
    [r_real,r_imag] = raiz.as_real_imag()
    grafROC.plot(r_real,r_imag,'o',color='orange',
             label ='polo:'+str(float(raiz))+';'+str(Q_raiz[raiz]))
grafROC.plot(1,0,'o',color='green',
         label ='radio:'+str(1))

plt.axhline(0,color='grey')
plt.axvline(0,color='grey')
plt.grid()
plt.legend()
plt.xlabel('Re[z]')
plt.ylabel('Imag[z]')
untitulo = 'ROC H[z]='+str(Xz)
plt.title(untitulo)

# h[n] usando tabla de transformadas
m = 10
n = sym.Symbol('n')
hn = 0.5833*((0.7)**n)*sym.Heaviside(n)+0.4166*((-0.5)**n)*sym.Heaviside(n)
fn = sym.lambdify(n,hn)
ki = np.arange(0,m,1)
xnk = fn(ki)

# grafica h[n]
figura, grafxn = plt.subplots()
plt.stem(ki,xnk)
plt.xlabel('ki')
plt.ylabel('h[n]')
plt.title('h[n]= '+str(hn))
plt.show()

s1Eva2009TII_T2 LTI DT bloques, respuesta impulso y paso

Ejercicio: 1Eva2009TII_T2 LTI DT bloques, respuesta impulso y paso

literal a. Respuesta al impulso, o H(s) del diagrama de bloques

a.1 H[z] Función de transferencia

Referencia: LTI DT – H[z] Diagrama de bloques con «1/z»

Dado que hay dos bloques (1/z) o de retraso , se identifica que los polinomios son de segundo orden (N=2). Los coeficientes de retraso se leen de arriba hacia abajo como a0,-a1,-a2 ordenados en la siguiente tabla:

coeficientes
retraso adelanto
a0 = 1 b0 = 1
a1 = -(3/4) = -3/4 b1 = 0
a2 = -(-1/8) = 1/8 b2 = 0

Se aplica lo mismo para los coeficientes de adelanto y se completa el modelo para los polinomios del numerador y denominador:

H[z] = \frac{b_0 z^N +b_1 z^{N-1} + \text{ ... } + b_{N-1} z + b_N}{z^N + a_1 z^{N-1} +\text{ ... } + a_{N-1}z +a_N} H[z] = \frac{z^2}{z^2 - \frac{3}{4} z+\frac{1}{8}}

a.2 H[z] en fracciones parciales

fracciones parciales modificadas,

\frac{H[z]}{z} = \Big( \frac{1}{z} \Big) \frac{z^2}{z^2 - \frac{3}{4} z+\frac{1}{8}} = \frac{z}{z^2 - \frac{3}{4} z+\frac{1}{8}}

El polinomio del denominador es:

Q[z] = z^2 - \frac{3}{4} z+\frac{1}{8}

con raíces en 1/4 y 1/2

con lo que se obtiene la expresión para usar el método de «cubrir» de Heaviside

\frac{H[z]}{z} = \frac{z}{(z-1/4)(z-1/2)}= \frac{k_1}{z-1/4}+ \frac{k_2}{z-1/2} k_1 = \frac{z}{\cancel{(z-1/4)}(z-1/2)}\Big|_{z=1/4} = \frac{1/4}{(1/4-1/2)} = -1 k_2 = \frac{z}{(z-1/4)\cancel{(z-1/2)}}\Big|_{z=1/2} = \frac{1/2}{(1/2-1/4)} = 2 \frac{H[z]}{z} = \frac{-1}{z-1/4}+ \frac{2}{z-1/2}

fracciones parciales

H[z] = \frac{-z}{z-1/4}+ \frac{2z}{z-1/2}

Revisando el resultado con el algoritmo en Python:

 Hz:
          2        
         z         
-------------------
 2                 
z  - 0.75*z + 0.125

 Hz/z:
           z           
-----------------------
     2                 
1.0*z  - 0.75*z + 0.125

 Hz/z.apart:
       1             2.0    
- ------------ + -----------
  1.0*z - 0.25   1.0*z - 0.5

 Hz = (Hz/z)*z
       z            2.0*z   
- ------------ + -----------
  1.0*z - 0.25   1.0*z - 0.5

 polos:     {0.250000000000000: 1, 0.500000000000000: 1}

 polos Re:  [1/4, 1/2]

Se obtiene h[n] usando la tabla de transformadas z

 Hz = (Hz/z)*z
       z            2.0*z   
- ------------ + -----------
  1.0*z - 0.25   1.0*z - 0.5
h[n] = - \Big(\frac{1}{4}\Big)^n \mu[n]+2\Big( \frac{1}{2} \Big)^n \mu[n]

a.3. Estabilidad del sistema

Los polos son menores que un radio de 1.
La ROC para el primer polo |z|>1/4 y para el segundo polo es |z|>1/2

Observaciones:

Las raíces características o frecuencias naturales del sistema se encuentran dentro del círculo de radio unitario. El sistema es asintóticamente estable, que implica que es BIBO estable.

h[n] no es de la forma k δ[n], por lo que el sistema global es con memoria.

a.4 Instrucciones con Python

Usando los algoritmo de la sección LTID Transformada z – X[z] Fracciones parciales modificadas con Python se adapta para el ejercicio

# Transformada z- Fracciones parciales
# Polos únicos, repetidos y complejos
# Lathi Ejercicio 5.3c pd495
# blog.espol.edu.ec/telg1001
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

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

Pz = z**2
Qz = z**2-(3/4)*z+1/8

# PROCEDIMIENTO
P = Pz.as_poly(z)
Q = Qz.as_poly(z)
Xz = Pz/Qz

# Raices Denominador
Q_raiz = sym.roots(Q)
Q_raizRe = sym.real_roots(Q)
rcompleja = len(Q_raiz)-len(Q_raizRe)

# Raices reales, únicas y repetidas
if (rcompleja<=0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z

def parametro_cuadratico(untermino):
    unparametro ={}
    # revisa denominador cuadratico
    [numerador,denominador] = (untermino).as_numer_denom()
    gradoD = 0
    coeficientesD = denominador
    gradoN = 0
    coeficientesN = numerador
    if not(denominador.is_constant()):
        denominador = denominador.as_poly()
        gradoD = denominador.degree()
        coeficientesD = denominador.coeffs()
    if not(numerador.is_constant()):
        numerador = numerador.as_poly()
        gradoN = numerador.degree()
        coeficientesN = numerador.coeffs()
    if gradoD == 2 and gradoN==2:
        a = float(coeficientesD[1])/2
        gamma2 = float(coeficientesD[2])
        gamma = np.sqrt(gamma2)
        A = float(coeficientesN[0])
        B = float(coeficientesN[1])
        rN = (A**2)*gamma2 + B**2 - 2*A*a*B
        rD = gamma2 - a**2
        r = np.sqrt(rN/rD)
        beta = np.arccos(-a/gamma)
        thetaN = A*a-B
        thetaD = A*np.sqrt(gamma2-a**2)
        theta = np.arctan(thetaN/thetaD)
        unparametro = {'r':r,
                       'gamma':gamma,
                       'beta':beta,
                       'theta':theta}
    return (unparametro)

# Terminos cuadraticos
parametros = [] # denominador cuadratico
if (rcompleja>0 and len(Q_raizRe)>0):
    # fracciones parciales modificadas
    Xzz = (P/Q)/z
    Xzm = Xzz.apart()
    # fracciones parciales restaurada
    terminos = Xzm.args
    Xzp = 0*z
    for untermino in terminos:
        Xzp = Xzp + untermino*z
        unparam = parametro_cuadratico(untermino*z)
        if len(unparam)>0:
            parametros.append(unparam)
if (rcompleja>0 and len(Q_raizRe)==0):
    Xzp = P/Q
    parametros.append(parametro_cuadratico(P/Q))

# SALIDA
print('\n Hz:')
sym.pprint(Xz)
if (len(Q_raizRe)>0):
    print('\n Hz/z:')
    sym.pprint(Xzz)
    print('\n Hz/z.apart:')
    sym.pprint(Xzm)
    print('\n Hz = (Hz/z)*z')
    sym.pprint(Xzp)
    print('\n polos:    ', Q_raiz)
print('\n polos Re: ', Q_raizRe)
if len(parametros)>0:
    for unparam in parametros:
        print('parametros cuadraticos: ')
        for cadauno in unparam.keys():
            print(cadauno,'\t',unparam[cadauno])


# grafica plano z imaginario
figura, grafROC = plt.subplots()
# limite
radio1 = plt.Circle((0,0),1,color='lightgreen',
                    fill=True)
radio2 = plt.Circle((0,0),1,linestyle='dashed',
                    color='lightgreen',fill=False)
grafROC.add_patch(radio1)
for raiz in Q_raiz.keys():
    [r_real,r_imag] = raiz.as_real_imag()
    radio_raiz = np.abs(raiz)
    nROC = plt.Circle((0,0),radio_raiz,
                      color='red',fill=True)
    grafROC.add_patch(nROC)
grafROC.add_patch(radio2) # borde r=1
grafROC.axis('equal')
# marcas de r=1 y valor a
for raiz in Q_raiz.keys():
    [r_real,r_imag] = raiz.as_real_imag()
    grafROC.plot(r_real,r_imag,'o',color='orange',
             label ='polo:'+str(float(raiz))+';'+str(Q_raiz[raiz]))
grafROC.plot(1,0,'o',color='green',
         label ='radio:'+str(1))

plt.axhline(0,color='grey')
plt.axvline(0,color='grey')
plt.grid()
plt.legend()
plt.xlabel('Re[z]')
plt.ylabel('Imag[z]')
untitulo = 'ROC H[z]='+str(Xz)
plt.title(untitulo)

# h[n] usando tabla de transformadas
m = 10
n = sym.Symbol('n')
hn = 0.5833*((0.7)**n)*sym.Heaviside(n)+0.4166*((-0.5)**n)*sym.Heaviside(n)
fn = sym.lambdify(n,hn)
ki = np.arange(0,m,1)
xnk = fn(ki)

# grafica h[n]
figura, grafxn = plt.subplots()
plt.stem(ki,xnk)
plt.xlabel('ki')
plt.ylabel('h[n]')
plt.title('h[n]= '+str(hn))
plt.show()

literal b. La respuesta de paso s[n]

La señal de entrada será x[n] = μ[n], con transformada z de:

X[z] = \frac{z}{z-1}

b.1 Y[z] ante entrada X[z] y con función de transferencia H[z]

Usando el resultado de la respuesta al impulso,

H[z] = \frac{-z}{z-1/4}+ \frac{2z}{z-1/2}

La salida Y[n] del sistema será:

Y[z] = X[z] H[z] = \frac{z}{z-1} \Bigg[ \frac{-z}{z-1/4} + \frac{2z}{z-1/2} \Bigg]

b.2 Y[z] en fracciones parciales

Fracciones parciales modificadas

\frac{Y[z]}{z} = \frac{1}{z-1} \Bigg[ \frac{-z}{z-1/4} + \frac{2z}{z-1/2} \Bigg] \frac{Y[z]}{z} = \frac{1}{(z-1)} \frac{-z}{z-1/4} + \frac{1}{(z-1)}\frac{2z}{z-1/2}

La expresión se puede operar por partes:

\frac{-z}{(z-1)(z-1/4)} = \frac{k_1}{z-1} + \frac{k_2}{z-1/4} k_1 = \frac{-z}{\cancel{(z-1)}(z-1/4)} \Big|_{z=1} = \frac{-1}{(1-1/4)}=-\frac{4}{3} k_2 = \frac{-z}{(z-1)\cancel{(z-1/4)}} \Big|_{z=1/4} = \frac{-1/4}{(1/4-1)} = \frac{1}{3}

y la segunda parte,

\frac{2z}{(z-1)(z-1/2)} = \frac{k_3}{z-1} + \frac{k_4}{z-1/2} k_3 = \frac{2z}{\cancel{(z-1)}(z-1/2)} \Big|_{z=1} = \frac{2(1)}{1-1/2} = 4 k_4 = \frac{2z}{(z-1)\cancel{(z-1/2)}} \Big|_{z=1/2} = \frac{2(1/2)}{1/2-1} = -2

Uniendo las partes de fracciones parciales modificadas se tiene:

\frac{Y[z]}{z} = \frac{-4/3}{z-1} + \frac{1/3}{z-1/4} + \frac{4}{z-1}+\frac{-2}{z-1/2}

restaurando las fraciones parciales

Y[z] = \frac{8}{3}\frac{z}{z-1} + \frac{1}{3}\frac{z}{z-1/4} -2\frac{z}{z-1/2}

Usando la tabla de transformadas z

y[n] = \frac{8}{3}\mu [n] + \frac{1}{3}\Big( \frac{1}{4}\Big)^n \mu [n] -2 \Big(\frac{1}{2}\Big) ^n \mu[n]

Usando el algoritmo en Python, se tiene:

 Yz:
              3             
             z              
----------------------------
(z - 1)*(z - 0.5)*(z - 0.25)

 Yz/z:
                 2                
                z                 
----------------------------------
     3         2                  
1.0*z  - 1.75*z  + 0.875*z - 0.125

 Yz/z.apart:
0.333333333333333       2.0       2.66666666666667
----------------- - ----------- + ----------------
   1.0*z - 0.25     1.0*z - 0.5        z - 1      

 Yz = (Yz/z)*z
0.333333333333333*z      2.0*z      2.66666666666667*z
------------------- - ----------- + ------------------
    1.0*z - 0.25      1.0*z - 0.5         z - 1       

 polos:  {0.2500: 1, 0.5000: 1, 1.0000: 1}

gráfica de salida y[n]

s1Eva2010TII_T3 LTI CT H(s) de bloques en paralelo y retraso

Referencia: 1Eva2010TII_T3 LTI CT H(s) de bloques en paralelo y retraso

El diagrama de bloques del sistema consta de dos sistema de primer orden en paralelo y el resultado en serie con un bloque retraso en tiempo.

La ecuación de respuesta a impulso H(s), siguiendo el diagrama se presenta como:

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

literal a

Para el desarrollo analítico se simplifica el problema, separando el bloque de atraso para el final, dado que el sistema es LTI los desplazamietos en el tiempo de la entrada tendrán semejante respuesta en la salida.

La función de transferencia H(s) o respuesta al impulso será:

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

Los componentes por ser de primer orden y estar en paralelo, no se requiere aplicar fracciones parciales.

Los polos del sistema se encuentran en el lado izquierdo del plano, por lo que sus componentes en el tiempo son decrecientes.

{polos,veces}:  {-2: 1, -3: 1}
 polos reales:  2  complejos:  0
 sobre lado derecho RHP: 0
 sobre Eje Imaginario, repetidos:  0  unicos: 0
 asintoticamente:  estable

NO existen polos en el lado derecho del plano, por lo que el sistema es asintoticamente estable, en consecuencia BIBO estable.

s1Eva2010TII_T3_polos H(s)literal b

Para la transformada inversa de Laplace, se recuerda que se tiene un componente de retraso en cascada, por lo que se ajusta H(s)*retraso en fracciones parciales, se presentan dos componentes:

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

que tienen la forma decreciente de la función, Para la gráfica sympy usa como θ(t)=μ(t)=Heaviside(t)

s1Eva2010TII_T3 x(t)h(t)=y(t)

literal c

A partir del diagrama de bloques se tiene:

H(s) =\frac{Y(s)}{X(s)}= \Big[ \frac{2s+5}{s^2+5s+6} \Big] e^{-2s} (s^2+5s+6)Y(s)=(2s+5)e^{-2s} X(s) s^2Y(s)+5sY(s)+6Y(s)=(2sX(s)+5X(s))e^{-2s}

que al convertir al dominio del tiempo se escribe como:

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

literal d

la respuesta ante la entrada escalon requiere que x(t) = μ(t), no se especifican condiciones iniciales del problema, por lo que se asumen iguales a cero. De la tabla de transformadas se obtiene X(s) = 1/s

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

que separando en fracciones parciales se convierte en:

Y(s) =\Big[\frac{5/6}{s} -\frac{1/3}{s+3}-\frac{1/2}{s+2} \Big] e^{-2s} y(t) =\Big[\frac{5}{6} -\frac{1}{3}e^{-3(t-2)}-\frac{1}{2}e^{-2(t-2)}\Big] \mu (t-2)

el resultado usando algoritmos es:

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

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

polosceros:
exp(-2*s) : {'Q_polos': {-2: 1, -3: 1}, 'P_ceros': {-5/2: 1}, 'Hs_k': (2*s + 5)/((s + 2)*(s + 3))}
Q_polos : {-2: 1, -3: 1}
P_ceros : {-5/2: 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

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

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

 Y(s)_total = ZIR + ZSR:
/      1           1        5 \  -2*s
|- --------- - --------- + ---|*e    
\  3*(s + 3)   2*(s + 2)   6*s/      

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

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) y estabilidad
Hs = (1/(s+2)+1/(s+3))*sym.exp(-2*s)
#Hs = 1+0*s cuando es constante

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

# 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 = 5
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(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) --------
muestras_H = 201
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,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()

s1Eva2010TII_T1 LTI CT Determinar h(t) con entrada x(t) y salida y(t)

Ejercicio: 1Eva2010TII_T1 LTI CT Determinar h(t) con entrada x(t) y salida y(t)

Dadas las señales de entrada y salida del sistema:

y(t) = \Big( 4 e^{-2t} + 6 e^{-3t} \Big) \mu (t) x(t) = 2 e^{-2t} \mu (t)

literal a. h(t)

Se usan las transformadas de Laplace para plantear la respuesta al impulso como  H(s) = Y(s)/X(s)

Y(s) = 4 \frac{1}{s+2} + 6 \frac{1}{s+3} = \frac{4(s+3)+6(s+2)}{(s+2)(s+3)} Y(s) = \frac{4s+12+6s+12}{(s+2)(s+3)} =\frac{10s+24}{(s+2)(s+3)} X(s) = 2 \frac{1}{s+2}

se obtiene la expresión para H(s),

H(s) = \frac{Y(s)}{H(s)} = \frac{\frac{10s+24}{(s+2)(s+3)}}{2 \frac{1}{s+2}} = \frac{2(5s+12)(s+2)}{2(s+2)(s+3)} H(s) = \frac{5s+12}{s+3}

aplicando fracciones parciales, conociendo que el grado M del polinomio P es igual al grado N del polinomio Q, existe un término constante mas terminos cero(k)/(s+polo(k)). El término constante es el coeficiente de s de mayor grado del numerador.

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

Aplicando la Transformada Inversa de Laplace, se tiene:

h(t) = 5\delta(t)-3e^{-3t} \mu (t)

con lo que comparando el resultado con la ecuación sugerida, a= 5, b=-3, c =3.

literal b. ecuación diferencial de h(t)

tomando la forma:

H(s) = \frac{5s+12}{s+3} \frac{Y(s)}{X(s)} = \frac{5s+12}{s+3} (s+3)Y(s) = (5s+12)X(s) sY(s)+3Y(s) = 5sX(s)+12X(s) \frac{\delta}{\delta t}y(t)+3y(t) = 5\frac{\delta}{\delta t}x(t)+12x(t)

Literal c. Diagrama de bloques

Literal d. polos

El sistema tiene un polo en -3 , que se encuentra en el lado izquierdo del plano s, lo que indica que tiene términos decrecientes y es asintóticamente estable. Lo mismo aplica para el caso de BIBO estable.

Tarea. literal e.

Para la señal indicada, considere revisar la propiedad de diferenciación en la tabla de propiedades transformada de Laplace.

\frac{\delta}{\delta t} x(t) \rightarrow sX(s) -x(0\text{–})

Resultados con Python

 y(t)
/   -2*t      -3*t\             
\4*e     + 6*e    /*Heaviside(t)
Y(s)
  6       4  
----- + -----
s + 3   s + 2

 x(t)
   -2*t             
2*e    *Heaviside(t)
X(s)
  2  
-----
s + 2

H(s) = P(s)/Q(s):
      3  
5 - -----
    s + 3

H(s) en factores:
5*s + 12
--------
 s + 3  

 h(t) :
                     -3*t             
5*DiracDelta(t) - 3*e    *Heaviside(t)

polosceros:
Q_polos : {-3: 1}
P_ceros : {-12/5: 1}

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

gráfica de H(s) con polos

s1Eva2010TII_T1_polos H(s)

gráfica de x(t),y(t), h(t)

s1Eva2010TII_T1_xh_y

Instrucciones en Python

# 1Eva2010TII_T1 LTI CT Determinar h(t) con entrada x(t) y salida y(t)
# Transformadas de Laplace, funciones
# http://blog.espol.edu.ec/telg1001/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1001 as fcnm

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

yt = (4*sym.exp(-2*t)+6*sym.exp(-3*t))*u

xt = 2*sym.exp(-2*t)*u

# Grafica, intervalo tiempo [t_a,t_b]
t_a = 0 ; t_b = 5
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Ys = fcnm.laplace_transform_suma(yt)
Xs = fcnm.laplace_transform_suma(xt)

# H(s) respuesta a estado cero
Hs = Ys/Xs
Hs = fcnm.apart_exp(Hs)
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(Q_polos)

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)

# SALIDA
print(' y(t)')
sym.pprint(yt)
print('Y(s)')
sym.pprint(Ys)
print('\n x(t)')
sym.pprint(xt)
print('X(s)')
sym.pprint(Xs)
print('\nH(s) = P(s)/Q(s):')
sym.pprint(Hs)
print('\nH(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])

# Graficas polos, H(s), con polos h(t) --------
muestras_H = 101
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,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()

s1Eva2009TII_T3 LTI CT y(t) desde h(t) y x(t) con términos escalón desplazados

Ejercicio: 1Eva2009TII_T3 LTI CT y(t) desde h(t) y x(t) con términos escalón desplazados

La función de transferencia h(t) y la señal de entrada x(t) son,

x(t) = 2 μ(t-1) – 2 μ(t-3)

h(t) = μ(t+1) – 2 μ(t-1) + μ(t-3)

La primera observación es que h(t) h(t) tiene un componente No causal, que se adelanta en tiempo μ(t+1) a x(t). La transformada unilateral de Laplace no aplica para este término, por lo que para el algoritmo se usa la entrada X(s) y se evita el error en la transformada. Por lo demás el algoritmo funciona bien.
Tarea: Revisar y justificar

Usando la función de transferencia h(t) con las transformadas de Laplace H(s)

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

simplificando la expresión como F(s)*(terminos exponenciales)

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

tiene la forma gráfica,

La señal de entrada usando las transformadas de Laplace X(s)

X(s) = 2\frac{1}{s}e^{-s} - 2\frac{1}{s}e^{-3s} X(s) = \frac{2}{s} \Big[e^{-s} - e^{-3s} \Big]

La señal de salida Y(s)=H(s)*X(s)

Y(s) = \frac{2}{s^2} \Big[ e^{s} - 2e^{-s}+ e^{-3s} \Big] \Big[e^{-s} - e^{-3s} \Big]

multiplicando los términos de exponenciales

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

usando la tabla de transformadas de Laplace en forma inversa:

y(t) = 2 (t) \mu (t) -6 (t-2) \mu (t-2) +6(t-4) \mu (t-4) -2(t-6) \mu(t-6) \Big]

s1Eva2009TII_T3_LTI xh_y

La respuesta al impulso H(s)no tiene polos en el lado derecho del plano s, por lo que las salidas son acotadas, en consecuencia el sistema es asintoticamente estable y también BiBO estable.

s1Eva2009TII_T3_LTI polos H(s)


usando el algoritmo, considerando que las condiciones iniciales son cero,:

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

 h(t) :
Heaviside(t - 3) - 2*Heaviside(t - 1) + Heaviside(t + 1)

polosceros:
exp(s) : {'Q_polos': {0: 1}, 'P_ceros': {}, 'Hs_k': 1/s}
exp(-3*s) : {'Q_polos': {0: 1}, 'P_ceros': {}, 'Hs_k': 1/s}
exp(-s) : {'Q_polos': {0: 1}, 'P_ceros': {}, 'Hs_k': -2/s}
Q_polos : {0: 1}
P_ceros : {}

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

 X(s): 
   -s      -3*s
2*e     2*e    
----- - -------
  s        s   

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

 ZSR respuesta estado cero:
ZSR :
        -2*s      -4*s      -6*s
2    6*e       6*e       2*e    
-- - ------- + ------- - -------
 2       2         2         2  
s       s         s         s   
yt_ZSR :
2*t*Heaviside(t) + (12 - 6*t)*Heaviside(t - 2) + (12 - 2*t)*Heaviside(t - 6) +
 (6*t - 24)*Heaviside(t - 4)

 Y(s)_total = ZIR + ZSR:
        -2*s      -4*s      -6*s
2    6*e       6*e       2*e    
-- - ------- + ------- - -------
 2       2         2         2  
s       s         s         s   

 y(t)_total = ZIR + ZSR:
2*t*Heaviside(t) + (12 - 6*t)*Heaviside(t - 2) + (12 - 2*t)*Heaviside(t - 6) +
 (6*t - 24)*Heaviside(t - 4)
>>>

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
Hs = sym.exp(s)/s - 2*sym.exp(-s)/s + sym.exp(-3*s)/s

# X(s) Señal de entrada
xt = 2*u.subs(t,t-1) - 2*u.subs(t,t-3)

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

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

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(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) --------
muestras_H = 101
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,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()

s1Eva2009TII_T1 LTI CT diagrama canónico y respuesta a impulso

Ejercicio1Eva2009TII_T1 LTI CT diagrama canónico y respuesta a impulso

a. A partir del diagrama canónico se puede reordenar, identificando las partes del numerador P(s) y denominador Q(s) para escribir la función de transferencia;

reordenando diagrama:

escribiendo la ecuación como:

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

Se usan los polos de la expresión para el denominador Q, para realizar las fracciones parciales:

 {polos,veces}:  {-1: 1, -2: 1}
H(s) = \Big[ \frac{3s+2}{(s+2)(s+1)} \Big] e^{-3s}

Al enfocarse en H(s) sin usar el retraso e-3s se tiene H1(s)

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

y se usa el método de «cubrir» de Heaviside para encontrar las constantes,

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

reemplazando las constantes en H1(s) o mejor en H(s) se obtiene expresiones mas simples en fracciones parciales.

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

b. Usar la tabla de transformadas de Laplace, y la propiedad de desplazamiento en t

h(t) = 4e^{-2(t-3)}\mu(t-3)- e^{-(t-3)}\mu(t-3)

Tarea: usar la forma analítica para encontrar la respuesta.

la gráfica de respuesta del sistema ht(t) y la salida y(t) del sistema ante una entrada impulso, serán iguales:

s1Eva2009TII_T1 xh_y

c. ¿Qué puede decir acerca de la estabilidad interna y externa?

Los polos se encuentran en el lado izquierdo del plano, sus componentes son exponenciales decrecientes. Por lo que se considera asintóticamente estable, tambien es acotado o BIBO estable.

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

s1Eva2009TII_T1_polosrevisando el resultado de los polos en conjunto con H(s)

s1Eva2009TII_T1_polos_ht

Resultados con el algoritmo:

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

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

polosceros:
exp(-3*s) : {'Q_polos': {-1: 1, -2: 1}, 
             'P_ceros': {-2/3: 1}}
Q_polos : {-1: 1, -2: 1}
P_ceros : {-2/3: 1}

Estabilidad de H(s):
 n_polos_real : 2
 n_polos_imag : 0
 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 :
/  4       1  \  -3*s
|----- - -----|*e    
\s + 2   s + 1/      
yt_ZSR :
/   3  -t      6  -2*t\                 
\- e *e   + 4*e *e    /*Heaviside(t - 3)

 Y(s)_total = ZIR + ZSR:
/  4       1  \  -3*s
|----- - -----|*e    
\s + 2   s + 1/      

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

Instrucciones en Python

Para resolver el ejercicio, se define la funcion H(s) por los polinomios P y Q, la entrada se define como un impulso δ(t).

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) y estabilidad
Hs = (3*s+2)/(s**2+3*s+2) * sym.exp(-3*s)
#Hs = 2*(s+5)/((s+4)*(s+3)) * sym.exp(-3*s)

# X(s) Señal de entrada
xt = d #+ d.subs(t,t-2)

# 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 = 0 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(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) --------
muestras_H = 101
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,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()

Referencia: LTI CT Laplace – H(s) Estabilidad del sistema con Python
Ejemplo 3, LTI CT Laplace – H(s) Diagramas de bloques y ecuaciones diferenciales

s1Eva2016TII_T4 LTI DT Sistema 2do orden

Ejercicio: 1Eva2016TII_T4 LTI DT Sistema 2do orden

literal a. Ecuación de diferencias de coeficientes constantes

Al diagrama del ejercicio se añaden las referencias de y[n], y[n-1] y y[n-2] de acuerdo a los bloques de retraso (delay) para escribir la expresión:

y[n] = x[n] + \frac{3}{4} y[n-1] - \frac{1}{8} y[n-2] y[n] - \frac{3}{4} y[n-1] + \frac{1}{8} y[n-2] = x[n]

Para usar los operadores ‘E’ y encontrar el polinomio característico, dado que el sistema es LTI – DT, se  puede desplazar 2 unidades:

y[n+2] - \frac{3}{4} y[n+1] + \frac{1}{8} y[n] = x[n+2]

La expresión usando notación de operadores ‘E’ es:

\Big[ E^2 - \frac{3}{4} E + \frac{1}{8} \Big] y[n] = E^2 x[n]

donde el numerador P[E] = E2
y el denominador Q[E] = E2 – (3/4)E +(1/8)

H(E) = \frac{P[E]}{Q[E]} = \frac{E^2}{E^2-\frac{3}{4}E+\frac{1}{8}}

literal b. Respuesta a impulso

A partir de la notación de operadores’E’,  se busca los modos característicos en el polinomio del denominador Q(E). También se expresa como y[n] como resultado de entrada cero x[n]=0

Q[E] = E^2 - \frac{3}{4} E + \frac{1}{8} \gamma^2 - \frac{3}{4} \gamma + \frac{1}{8} = 0

usando la fórmula o el algoritmo se busca las raices del polinomio Q(E):

(\gamma - \frac{1}{4})(\gamma - \frac{1}{2}) = 0 \gamma_1 = \frac{1}{4} ; \gamma_2 = \frac{1}{2}

Siendo raices reales y no repetidas, la respuesta del sistema tiene la forma:

y_c[n] = c_1 \Big( \frac{1}{4} \Big) ^n +c_2 \Big( \frac{1}{2} \Big) ^n

se convierte a la forma:

h[n] =\frac{b_n}{a_n} \delta [n] + y_c[n] \mu [n] h[n] = \Bigg[ c_1 \Big( \frac{1}{4} \Big) ^n +c_2 \Big( \frac{1}{2} \Big) ^n \Bigg] \mu [n]

Para encontrar los coeficientes c1 y c2, se requieren dos valores iniciales. En el literal a se indica que el sistema es causal, en consecuencia «No es posible obtener una salida antes que se aplique la entrada» y tenemos que: y[-1] = 0; y[-2] = 0.

La ecuación original de y[n] para respuesta a impulso tiene entrada x[n]=δ[t]

y[n] = \frac{3}{4} y[n-1] - \frac{1}{8} y[n-2] + x[n] y[n] = \frac{3}{4} y[n-1] - \frac{1}{8} y[n-2] + \delta [n]

donde el impulso tiene valor de 1 solo en n=0, haciendo equivalente h[n] equivale a y[n] , entonces:

h[n] = \delta [n] + \frac{3}{4} h[n-1] - \frac{1}{8} h[n-2]

n = 0

h[0] = \delta [0] + \frac{3}{4} h[-1] - \frac{1}{8} h[-2] h[0] = \delta [0] + \frac{3}{4} (0) - \frac{1}{8} (0) = 1

n=1

h[1] = \delta [1] + \frac{3}{4} h[0] - \frac{1}{8} h[-1] h[1] = (0) + \frac{3}{4} (1) - \frac{1}{8} (0) = \frac{3}{4}

con lo que es posible encontrar c1 y c2,

h[0] = 1 = c_1 \Big( \frac{1}{4} \Big) ^0 \mu [0] + c_2 \Big( \frac{1}{2} \Big) ^0 \mu (0) c_1 + c_2 = 1 h[1] = \frac{3}{4} = c_1 \Big( \frac{1}{4} \Big) ^1 \mu [1] + c_2 \Big( \frac{1}{2} \Big) ^1 \mu (1) \frac{1}{4} c_1 + \frac{1}{2} c_2 = \frac{3}{4}

resolviendo:

c_1 = -1 c_2 = 2 h[n] = - \Big( \frac{1}{4} \Big) ^n \mu [n] + 2 \Big( \frac{1}{2} \Big) ^n \mu [n] h[n] = \Bigg[ 2 \Big( \frac{1}{2} \Big) ^n - \Big( \frac{1}{4} \Big) ^n \Bigg] \mu [n]

Observaciones:

Las raíces características o frecuencias naturales del sistema se encuentran dentro del círculo de radio unitario. El sistema es asintóticamente estable, que implica que es BIBO estable.

h[n] no es de la forma k δ[n], por lo que el sistema global es con memoria.

La forma de respuesta al impulso se vuelve evidente que el sistema es IIR.

Algoritmo en Python

Usando el algoritmo de LTID – Respuesta impulso. Ejercicio con Python:

 respuesta impulso: 
hi: [1.   0.75]
Matriz: 
[[1.   1.  ]
 [0.5  0.25]]
Ch:  [ 2. -1.]
n>=0, hn:  
      n          n
- 0.25  + 2.0*0.5 
>>> 

Instrucciones Python

# Sistema LTID. Respuesta impulso
# QE con raices Reales NO repetidas Lathi ejemplo 3.13 pdf271
# Lathi ejemplo 3.18 y 3.19 pdf278,
# blog.espol.edu.ec/telg1001
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
# coeficientes E con grado descendente
QE = [1.0, -3/4, 1/8]
PE = [1., 0., 0.]
# condiciones iniciales ascendente ...,y[-2],y[-1]
inicial = [0, 0.]

tolera = 1e-6  # casi_cero
muestras = 10  # para grafica

# PROCEDIMIENTO
# Respuesta a ENTRADA CERO
# raices, revisa numeros complejos
gamma = np.roots(QE)
revisaImag = np.iscomplex(gamma)
escomplejo = np.sum(revisaImag)

# coeficientes de ecuacion
m_q = len(QE)-1
Ac = np.zeros(shape=(m_q,m_q),dtype=float)

# revisa si parte compleja <tolera o casi_cero 
if escomplejo>0:
    for i in range(0,m_q,1):
        valorimag = np.imag(gamma[i])
        if np.abs(valorimag)<tolera:
            gamma[i] = float(np.real(gamma[i]))
    sumaj = np.sum(np.abs(np.imag(gamma)))
    if sumaj <tolera:
        print(sumaj)
        gamma = np.real(gamma)
        escomplejo = 0

# revisa repetidos
unicoscuenta = np.unique(gamma,return_counts=True)
repetidas = np.sum(unicoscuenta[1]-1)

# Respuesta impulso h[n]
ki = np.arange(-m_q,m_q,1)
hi = np.zeros(m_q, dtype=float)
xi = np.zeros(2*m_q, dtype=float)
xi[m_q] = 1 # impulso en n=0
Ah = np.zeros(shape=(m_q,m_q),dtype=float)

# h[n] iterativo
p_n = len(PE)
for i in range(0,m_q,1):
    for k in range(m_q,2*m_q,1):
        hi[i] = hi[i] - QE[k-m_q]*hi[m_q-k+1]
    for k in range(0,p_n,1):
        hi[i] = hi[i] + PE[k]*xi[m_q+i]
Bh = np.copy(hi[:m_q])

# coeficientes de ecuacion
for i in range(0,m_q,1):
    for j in range(0,m_q,1):
        Ah[i,j] = gamma[j]**(i)
ch = np.linalg.solve(Ah,Bh)

# ecuacion hn
n = sym.Symbol('n')
hn = 0*n
for i in range(0,m_q,1):
    hn = hn + ch[i]*(gamma[i]**n)

# SALIDA
if escomplejo == 0: 
    print('\n respuesta impulso: ')
    print('Bh:',Bh)
    print('Matriz: ')
    print(Ah)
    print('Ch: ',ch)
    print('n>=0, hn:  ')
    sym.pprint(hn)
else:
    print(' existen raices con números complejos.')
    print(' usar algoritmo de la sección correspondiente.')
    

# grafica datos
ki = np.arange(-m_q,muestras,1)
hi = np.zeros(muestras+m_q)

if escomplejo == 0: 
    # evaluación de h[n]
    h_n = sym.lambdify(n,hn)
    hi[m_q:] = h_n(ki[m_q:])

    # grafica h[n]
    plt.stem(ki,hi,label='h[n]',
             markerfmt ='C1o',
             linefmt='C2--')
    plt.legend()
    plt.grid()
    plt.ylabel('h[n]')
    plt.xlabel('ki')
    plt.title('h[n] = '+str(hn))
    plt.show()


literal c.  respuesta de paso s[n]

s[n] = \sum_{k=-\infty}^{n} h[k] s[n \rightarrow \infty] = \sum_{k=-\infty}^{\infty} h[k] = \sum_{k=-\infty}^{\infty} \Bigg[ 2 \Big( \frac{1}{2} \Big) ^k - \Big( \frac{1}{4} \Big) ^k \Bigg] \mu [k] = \sum_{k=0}^{\infty} \Bigg[ 2 \Big( \frac{1}{2} \Big) ^k - \Big( \frac{1}{4} \Big) ^k \Bigg]

considerando que:

\sum_{k=0}^{\infty} \alpha ^k = \frac{1}{1-\alpha} ; |\alpha|<1

se reemplaza con:

s[n \rightarrow \infty] = \sum_{k=0}^{\infty} \Bigg[ 2 \Big( \frac{1}{2} \Big) ^k - \Big( \frac{1}{4} \Big) ^k \Bigg] = 2\frac{1}{1-\frac{1}{2}} - \frac{1}{1-\frac{1}{4}} = 2(2) - \frac{4}{3} s[n \rightarrow \infty] =\frac{8}{3}

un método mas detallado es:

s[n] = \sum_{k=-\infty}^{n} h[k] s[n] = \sum_{k=-\infty}^{n} \Bigg[ 2 \Big( \frac{1}{2} \Big) ^k - \Big( \frac{1}{4} \Big) ^k \Bigg] \mu [k] s[n] = \sum_{k=0}^{n} \Bigg[ 2 \Big( \frac{1}{2} \Big) ^k - \Big( \frac{1}{4} \Big) ^k \Bigg] s[n] = 2 \sum_{k=0}^{n} \Big( \frac{1}{2} \Big) ^k - \sum_{k=0}^{n} \Big( \frac{1}{4} \Big) ^k

considerando que:

\sum_{k=0}^{n} \alpha ^k = \frac{1-\alpha ^{n+1}}{1-\alpha}

se tiene lo siguiente,

s[n] = 2 \frac{1-\frac{1}{2}^{n+1}}{1-\frac{1}{2}} - \frac{1-\frac{1}{4}^{n+1}}{1-\frac{1}{4}} s[n] = 2 \frac{1-\frac{1}{2} \Big(\frac{1}{2} \Big)^{n}}{\frac{1}{2}} - \frac{1-\frac{1}{4} \Big( \frac{1}{4} \Big)^{n}}{\frac{3}{4}} s[n] = 4 \Bigg[ 1-\frac{1}{2} \Big(\frac{1}{2}\Big)^{n} \Bigg]- \frac{4}{3} \Bigg[ 1- \frac{1}{4} \Big(\frac{1}{4} \Big)^{n} \Bigg] s[n] = 4-2 \Big(\frac{1}{2}\Big)^{n} -\frac{4}{3}+ \frac{1}{3} \Big(\frac{1}{4} \Big)^{n} lim _ {n \to \infty}s[n] = lim _ {n \to \infty} \Bigg[ \frac{8}{3}-2 \Big(\frac{1}{2}\Big)^{n} + \frac{1}{3} \Big(\frac{1}{4} \Big)^{n} \Bigg] lim _ {n \to \infty}s[n] = \frac{8}{3}

que es el mismo resultado que con el método presentado al inicio del literal.



Solución alterna: Usando transformada z

Revisar: LTID Transformada z – X[z] Fracciones parciales modificadas con Python

A partir de la ecuación de diferencias:

y[n+2] - \frac{3}{4} y[n+1] + \frac{1}{8} y[n] = x[n+2]

en transformada z (unidad 7), que es semejante a operador ‘E’,

z^2 Y[z]- \frac{3}{4} z Y[z] + \frac{1}{8}Y[z] = z^2 X[z] \Big[ z^2 - \frac{3}{4} z + \frac{1}{8}\Big] Y[z] = z^2 X[z]

la función de transferencia o respuesta al impulso se expresa como:

H[z] = \frac{X[z]}{Y[z]} = \frac{z^2}{z^2 - \frac{3}{4} z + \frac{1}{8}}

usandolas raices del polinomio Q(E)

H[z] = \frac{z^2}{(z - \frac{1}{4})(z - \frac{1}{2})}

para facilitar las operaciones se usan fracciones parciales modificadas,

\frac{H[z]}{z} = \frac{z}{(z - \frac{1}{4})(z - \frac{1}{2})} = \frac{k_1}{z - \frac{1}{4}} + \frac{k_2}{z - \frac{1}{2}}

Usando el método de «cubrir» de Heaviside:

k_1 = \frac{z}{\cancel{(z - \frac{1}{4})}(z - \frac{1}{2})} \Big|_{z=1/4} = \frac{1/4}{1/4 - \frac{1}{2}} = -1 k_2 = \frac{z}{(z - \frac{1}{4})\cancel{(z - \frac{1}{2})}} \Big|_{z=1/2} = \frac{1/2}{1/2 - \frac{1}{4}} = 2

reemplazando k1, k2 y restaurando a fracciones parciales al multiplicar ambos lados por z,

H[z] = \frac{-z}{z - \frac{1}{4}} + \frac{2z}{z - \frac{1}{2}}

se puede usar la tabla de transformadas z para obtener el mismo resultado que el desarrolo anterior con operador ‘E’

h[n] = - \Big( \frac{1}{4} \Big) ^n \mu [n] + 2 \Big( \frac{1}{2} \Big) ^n \mu [n]

s1Eva2016TII_T3 LTI DT Sistema 1er orden con constante a

Ejercicio: 1Eva2016TII_T3 LTI DT Sistema 1er orden con constante a

A partir del diagrama se tiene que:

literal a.1 ecuación que relaciona entrada-salida

El sistema es discreto de primer orden,

y[n] = x[n] + a y[n-1] y[n] - a y[n-1] = x[n]

literal a.2 Respuesta al impulso

Para la respuesta de impulso se escribe la expresión en notación de operador E, con x[n]=δ[n]

y[n+1] - a y[n] = x[n+1] Ey[n] - a y[n] = E x[n] (E -a) y[n] = E x[n]

para determinan las raíces del Q[E]

Q[E] = (\gamma-a) = 0 \gamma = a

que establece la ecuación

y_c[n] = c_1 a^n

Para encontrar el valor de la constante se evalua la expresión del sistema con  n=0, teniendo en cuenta que la ecuación general de respuesta a impulso:

h[n] = \frac{b_N}{a_N} \delta (n) + y_c[n] \mu [n] h[n] = y_c[n] \mu [n] = c_1 a^n \mu [n] y[n] = x[n] +a y[n-1] h[n] = \delta [n] +a h[n-1] h[0] = \delta [0] +a h[0-1]

El sistema es causal si para n<0 los valores son cero, h[-1] =0

h[0] = 1 +a (0) = 1

con lo que  par encontrar la constante c1,

h[0] = 1 = c_1 a^0 \mu [0] =c_1 (1)(1) = c_1

la constante c1 =1, quedando com respuesta a impulso:

h[n] = a^n \mu [n]

literal a.3 Respuesta de paso

s[n] = \sum_{k=-\infty}^{n} h[k] = \sum_{k=-\infty}^{n} a^k \mu [k]

Siendo:

\sum_{k=0}^{n} \alpha^k = \frac{1-\alpha^{n+1}}{1-\alpha} \mu [n]

se tiene que:

s[n] = \frac{1-a^{n+1}}{1-a} \mu [n]

Solución alterna, El mismo resultado se obtiene mediante:

s[n] = h[n] \circledast \mu [n]

usando la tabla de convolucion discreta:

= \sum_{k=- \infty}^{\infty} a^k \mu[k] \mu[n-k] = \sum_{k=0}^{n} a^k s[n] = \frac{1-a^{n+1}}{1-a} \mu [n]

literal b.

dado que h[n] no es de la forma k δ[n], el sistema es con memoria.

Si las raíces características, valores característicos o frecuencias naturales de un sistema, se encuentran dentro del círculo de radio unitario, el sistema es asintóticamente estable e implica que es BIBO estable.

h[n] = an μ[n] el sistema es de tipo IIR.


literal c. Sistema inverso

h[n] \circledast h_{inv} = \delta[n] y[n] = x[n] + a y[n-1] h[n] = \delta[n] + a h[n-1] h[n] - a h[n-1] = \delta[n] h[n] \circledast \delta [n] - a h[n] \circledast \delta [n-1] = \delta[n] h[n] \circledast (\delta [n] - a \delta [n-1]) = \delta[n] h_{inv}[n] = (\delta [n] - a \delta [n-1]) = \delta[n]

h[n] no es de la forma k δ[n], por lo que el sistema inverso es con memoria.

h[n]=0 para n<0, se entiende que el sistema inverso es Causal.

La respuesta impulso del sistema inverso es absolutamente sumable o convergente, por lo que el sistema es BIBO estable.

la forma de hinv[n] el sistema es de tipo FIR.


literal d. Diagrama de bloques del sistema inverso

El diagrama del sistema inverso sería:

 

s1Eva2016TII_T2 LTI CT bloques en paralelo-serie con Laplace

Ejercicio: 1Eva2016TII_T2 LTI CT bloques en paralelo-serie con Laplace1Eva2012TII_T4 LTI CT bloques en paralelo-serie con Laplace, 1Eva2011TII_T3 LTI CT H(s) desde expresión con operadores D

Ejercicio resuelto: [ transformada de Laplace ] [integral convolución]

..


literal a. Funcion H(s) global

en el diagrama se muestran dos sistemas LTIC de primer orden, ambos sistemas se encuentran multiplicados por una constante y en paralelo.

La función de transferencia de un sistema LTIC de primer orden es
H(s) = 1/(s-a)

Con lo que el sistema se puede reescribir como:

H(s) = 3\frac{1}{s+7} + 12 \frac{1}{s-4}

Ceros y polos de H(s)

H(s) = \frac{3(s-4)+12(s+7)}{(s+7)(s-4)} = \frac{3s-12+12s+84}{(s+7)(s-4)} = \frac{15s+72}{(s+7)(s-4)}

los ceros se toman del numerador:

15s+72 = 0 s=\frac{-72}{15} = -\frac{24}{5}

los polos se toman del denominador:

(s+7)(s-4) = 0 s=-7 ; s= 4

s1Eva2016TII_T2_polos H(s)

literal b. La respuesta impulso h(t)

Los polos tienen signo diferente, lo que indica que el intervalo donde se encuentran incluyen al eje imaginario en el intervalo, por lo que el sistema es BIBO inestable. Por tener un polo a derecha del plano es asintoticamente inestable, su respuesta tiene componentes que crecen en el tiempo.

Como la región de convergencia ROC no se encuentra a la derecha de todos los polos, se tiene concluye que el sistema NO es causal.

h(t) = \mathscr{L}^{-1} [H(s)] = \mathscr{L}^{-1} \Big[ 3\frac{1}{s+7} + 12 \frac{1}{s-4}\Big] = \mathscr{L}^{-1} \Big[ 3\frac{1}{s+7} \Big] + \mathscr{L}^{-1} \Big[ 12 \frac{1}{s-4}\Big]

usando la tabla de transformadas o Sympy de Python

h(t) = 3e^{-7t} \mu (t) + 12 e^{4t} \mu (t)

s1Eva2016TII_T2_ht

literal c. Salida ante entrada x(t)

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

1Eva2016TII_T2 ZSR animado

X(s) = \frac{1}{s+5} Y(s) = X(s) H(s) =\frac{1}{s+5} \Bigg[\frac{15s+72}{(s+7)(s-4)}\Bigg] =\frac{15s+72}{(s+5)(s+7)(s-4)}

se usa fracciones parciales

=-\frac{3}{2(s+7)} + \frac{1}{6(s+5)}+\frac{4}{3(s-4)} y(t) = \mathscr{L}^{-1} [Y(s)] y(t) = \mathscr{L}^{-1} \Bigg[-\frac{3}{2(s+7)} + \frac{1}{6(s+5)}+\frac{4}{3(s-4)} \Bigg]

usando la tabla de transformadas de Laplace :

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

s1Eva2016TII_T2_polos y(t)


Algoritmo en Python

realizado a partir de LTIC Laplace – Algoritmo Python para analizar estabilidad H(s), Y(s) con entrada cero, estado cero, condiciones iniciales

 H(s) = P(s)/Q(s):
  3       12 
----- + -----
s + 7   s - 4
 H(s) en factores:
  3*(5*s + 24) 
---------------
(s - 4)*(s + 7)

 h(t) :
/    4*t      -7*t\             
\12*e    + 3*e    /*Heaviside(t)

polosceros:
Q_polos : {4: 1, -7: 1}
P_ceros : {-24/5: 1}

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

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

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

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

 Y(s)_total = ZIR + ZSR:
      3           1           4    
- --------- + --------- + ---------
  2*(s + 7)   6*(s + 5)   3*(s - 4)

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

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) y estabilidad
Hs = 3/(s+7) + 12/(s-4)
#Hs = 1+0*s cuando es constante

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

# 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 = -0.2 ; t_b = 0.5
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
Hs = fcnm.apart_exp(Hs) # fracciones parciales
Hs_fc = fcnm.factor_exp(Hs) # en factores
Hs_Qs2 = fcnm.Q_cuad_parametros(Hs_fc)

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

estable = fcnm.estabilidad_asintotica(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) --------
muestras_H = 201
figura_s  = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,f_nombre='H',solopolos=True)
figura_Hs = fcnm.graficar_Fs(Hs_fc,Q_polos,P_ceros,muestras=muestras_H,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()

Ejercicio resuelto: [ transformada de Laplace ] [integral convolución]
..


Ejercicio resuelto con integral de convolución, unidad 3

Algoritmo desarrollado con LTI CT Respuesta del Sistema Y(s)=ZIR+ZSR con Sympy-Python

 ZIR(t):
0

 h(t):
/    4*t      -7*t\             
\12*e    + 3*e    /*Heaviside(t)

 ZSR(t):
xh :
    4*t  -9*tau                                     
12*e   *e      *Heaviside(tau)*Heaviside(t - tau) + 

     -7*t  2*tau            
+ 3*e    *e     *Heaviside(tau)*Heaviside(t - tau)
xcausal : True
hcausal : True
[tau_a,tau_b] : [0, t]
intercambia : False
cond_graf : True
ZSR :
/   4*t    -5*t      -7*t\             
|4*e      e       3*e    |             
|------ + ----- - -------|*Heaviside(t)
\  3        6        2   /             

 y(t) = ZIR(t)+ZSR(t):
/   4*t    -5*t      -7*t\             
|4*e      e       3*e    |             
|------ + ----- - -------|*Heaviside(t)
\  3        6        2   /    

Las gráficas resultan iguales a las mostradas

Instrucciones con Python

# Respuesta total del sistema
# y(t) = ZIR(t) + ZSR(t)
# http://blog.espol.edu.ec/telg1001/lti-ct-yszirzsr-respuesta-del-sistema-con-sympy-python/
# Revisar causalidad de x(t) y h(t)
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import telg1001 as fcnm

# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) - 28*y(t)
RHS = 15*sym.diff(x(t),t,1,evaluate=False) + 72*x(t)
ecuacion = sym.Eq(LHS,RHS)

# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [0,0]

# entrada x(t)
x = sym.exp(-5*t)*u

# grafica intervalo [t_a,t_b]
t_a = 0; t_b = .5
muestras = 201

# PROCEDIMIENTO
# Respuesta entrada cero ZIR
sol_ZIR = fcnm.respuesta_ZIR(ecuacion,cond_inicio,t0)
ZIR = sol_ZIR['ZIR']

# Respuesta al impulso h(t)
sol_h = fcnm.respuesta_impulso_h(ecuacion)
h = sol_h['h']

# respuesta a estado cero ZSR
sol_ZSR = fcnm.respuesta_ZSR(x,h)
ZSR = sol_ZSR['ZSR']
xh  = sol_ZSR['xh']

# respuesta a y(t) = ZIR(t)+ZSR(t)
y_total = ZIR+ZSR

# revisa si grafica ZSR
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')
    
# SALIDA
print('\n ZIR(t):')
sym.pprint(ZIR)

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

print('\n ZSR(t):')
fcnm.print_resultado_dict(sol_ZSR)
if sol_ZSR['hcausal']==False:
    print('revisar causalidad de h(t)')
if sol_ZSR['xcausal']==False:
    print('revisar causalidad de x(t)')
if sol_ZSR['intercambia']:
    print('considere intercambiar h(t)con x(t)')
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')

print('\n y(t) = ZIR(t)+ZSR(t):')
sym.pprint(y_total)

# GRAFICA
figura_ZIR = fcnm.graficar_ft(ZIR,t_a,t_b,
                              muestras,'ZIR')
figura_h   = fcnm.graficar_ft(h,t_a,t_b,
                              muestras,'h')
# grafica animada de convolución
n_archivo = '' # sin crear archivo gif animado 
#n_archivo = '1Eva2016TII_T2_ZSR' # requiere 'imagemagick'
if sol_ZSR['cond_graf']:
    fig_ZSR = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
                                 muestras,y_nombre='ZSR')
    fig_ytotal = fcnm.graficar_xhy(ZIR,ZSR,y_total,t_a,t_b,
                                   muestras,x_nombre='ZIR',
                                   h_nombre='ZSR',y_nombre='y')
    figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,-t_b,t_b,
                      muestras, reprod_x = 4,y_nombre='ZSR',
                      archivo_nombre = n_archivo)
plt.show()