3.3.2 Respuesta a impulso – Diagrama Bloques

Para usar el simulador con diagrama de bloques para respuesta a impulso unitario, es necesario usar la expresión:

h(t) = b_0 \delta (t) +P(D)[y_n(t) \mu(t)]

Donde si el orden de P(D) es menor al orden de Q(D) entonces b0=0.

Para el ejercicio desarrollado en el sistema modelo, P(D)=D, por lo que hay una sola derivada, que se puede obtener de la entrada del primer integrador de y(t).

( D^2 + 3D +2 ) y(t) = Dx(t)

Para obtener h(t), el impulso unitario x(t), es semejante a tener una condicion inicial sin señal de entrada, siguiendo la teoría descrita de sobre condiciones iniciales. El impulso unitario representa que y'(t) inicia con 1, por ser el término de más alto orden en P(D) y se configura el diagrama de bloques:

Con lo que la respuesta del sistema al impulso unitario que se obtiene es:

El parámetros inicial se configuró en:

Como referencia para observación se muestra la salida y(t):

 

 

3.4.2 Respuesta a estado cero – Diagrama Bloques

Se continúa con el ejercicio propuesto en la sección Sistema-Modelo:

Referencia: Lathi ejemplo 2.6 pdf129

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

el sistema tiene respuesta a impulso:

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

El ejemplo es continuación del presentado para respuesta a entrada cero, que tiene la ecuación:

(D^2 + 3D +2)y(t) = Dx(t)

Al diagrama de Bloques desarrollado en la sección 3.2.2 Respuesta a entrada cero, se añaden los bloques para formar la señal de entrada. La variable tiempo se indica con el Bloque Reloj, el escalón μ(t) se configura para empezar en cero y los valores de -3 del exponente se los indica con un bloque de ganancia antes de eu.

Para diferenciación de los bloques, la entrada se usa el color azul, el sistema en color naranja y la salida en color morado.

Para continuar con la simulación, asegúrese que no existan condiciones iniciales en el sistema, es decir los valores de estado inicial para los bloques 1/s se encuentren en cero.

Los osciloscopios tiene su referencia de tiempo con otro reloj. El osciloscopio de la señal de entrada x(t) se muestra primero, y la salida luego.

la señal de respuesta se obtiene de forma semejante a la respuesta al impulso, es decir en el punto donde está la primera derivada. Revisar teoría para aplicar la forma: P(D)[y(t)μ(t)], que en éste ejercicio se aplica una D.

El resultado de y(t) solo como referencia es el que se obtiene en el osciloscopio ubicado a la derecha abajo.

3.3.1 Respuesta a impulso – Python

1. Fórmulas simbólicas con Sympy

Siguiendo las instrucciones de la sección teórica, usando los coeficientes de la ecuación y condiciones de inicio, se modifica el algoritmo para encontrar la respuesta a entrada cero. Las condiciones de inicio se establecen como y(0)=0 y y'(0)=0 para el problema planteado.

La ecuación diferencial: 
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
Solución general: 
           -t       -2*t
y(t) = C1*e   + C2*e    

Ecuaciones de condiciones iniciales
0 = C1 + C2
1 = -C1 - 2*C2
{C1: 1, C2: -1}
Solución particular:
        -t    -2*t
y(t) = e   - e    

Solución h(t): 
          -t      -2*t
h(t) = - e   + 2*e    
>>> 

se desarrolla con las instrucciones:

# Lathi 2.1.b pdf 116
# Respuesta a entrada cero con Sympy
# solucion para (D^2+ 3D + 2)y = 0
import sympy as sym

# INGRESO
# Condiciones iniciales
t0 = 0 ; y0 = 0 ; dy0 = 1

t = sym.Symbol('t')
y = sym.Function('y')
ecuacion = sym.Eq(sym.diff(y(t),t,t) + 3*sym.diff(y(t),t)+2*y(t),0)

# PROCEDIMIENTO
general = sym.dsolve(ecuacion, y(t))
general = general.expand()

# Condiciones iniciales
# rhs, righ hand side. Lado derecho de ecuacion
s0 = sym.Eq( y0, general.rhs.subs(t,t0))
s1 = sym.Eq(dy0, general.rhs.diff(t).subs(t,t0))
constante = sym.solve([s0,s1])

# Solucion particular
particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

h = sym.Function('h')
rimpulso = sym.Eq(h(t),particular.rhs.diff(t))

# SALIDA
print('La ecuación diferencial: ')
sym.pprint(ecuacion)
print('Solución general: ')
sym.pprint(general)

print('\nEcuaciones de condiciones iniciales')
sym.pprint(s0)               
sym.pprint(s1)
print(constante)
print('Solución particular:')
sym.pprint(particular)

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

y con gráfica:

# GRAFICA
import numpy as np
a = 0 ; b =5 ; muestras = 101
ht = sym.lambdify(t, rimpulso.rhs, 'numpy')

ti = np.linspace(a,b,muestras)
hi = ht(ti)

import matplotlib.pyplot as plt
plt.plot(ti,hi, color='red', label='y(t)')
plt.title('Respuesta impulso unitario: y(t)='+str(rimpulso.rhs))
plt.xlabel('t')
plt.grid()
plt.show()

2. Usando un sistema LTI definido con Scipy.signal

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación y condiciones de inicio, la respuesta a un impulso unitario δ(t) o función de transferencia h(t) se calcula con una versión muestreada xi a partir del tiempo muestreado en el intervalo [a,b].

Del problema del Sistema-Modelo y la ecuacion de segundo orden, se tiene dos opciones para la respuesta:

– scipy.signal.impulse() donde solo se da el vector ti
– aplicar la teoria de los coeficientes donde solo la derivada de mayor orden de P(D) evaluada en cero, tiene el valor de uno. Todos los demás términos tienen el valor de cero. Revisar la parte conceptual de la sección para respuesta a impulso.

Para el ejercicio y como comprobación se calculan las dos formas, con entrada cero y condicion inicial [1,0] y con la función de scipy, que son las usadas en la gráfica.

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)
import numpy as np
import scipy.signal as senal

# INGRESO
# tiempo [a,b]
a = 0 ; b = 5 ; muestras=101

# Señal de entrada
u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: np.piecewise(t,t==0,[1,0])
# 10*np.exp(-3*t)*u(t)

# Sistema
# coeficients Q  P de la ecuación diferencial
Q = [1., 3, 2]
P = [1., 0]
# condiciones de inicio [y'0,y0]
inicio = [1,0]

# PROCEDIMIENTO
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)
[t_y, yi, yc] = senal.lsim(sistema, 0,ti, inicio)
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.subplot(211)
plt.suptitle('Respuesta a Impulso Unitario')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')

plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(t_y, yi, color='magenta', label='yi(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='yi(t)')
# plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

3. Algoritmo de Runge-Kutta d2y/dx2 de Análisis numérico

Para el caso del algoritmo de Runge-Kutta, se aplica las condiciones iniciales de y(0)= 0 y y'(0)=1, siguiendo lo descrito en la parte teórica y se toman los valores de z=y’ para obtener el resultado.

Las instrucciones aplicadas al algoritmo de entrada cero se muestran a continuación:

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA

f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = 1


h = 0.1
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)

ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(ti,zi, color = 'Red', label='dy/dt')

plt.ylabel('dy/dt')
plt.xlabel('t')
plt.title('Respuesta impulso: Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Recuerde que se aplican las modificaciones de acuerdo al planteamiento del problema, por lo que revise la teoría antes de aplicar el algoritmo.

3.4.1 Respuesta a estado cero – Python

Se continúa con el ejercicio propuesto en la sección Sistema-Modelo:

Referencia:  Lathi ejemplo 2.6 pdf129

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

el sistema tiene respuesta a impulso:

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

El ejemplo es continuación del presentado para respuesta a entrada cero, que tiene la ecuación:

(D^2 + 3D +2)y(t) = Dx(t)

Se desarrollan el ejercicio de varias formas para obtener la gráfica de salida usando:

  1. Un sistema con Scipy.signal.lti() y simulado con Scipy.signal.lsim()
  2. La integral de convolución con Numpy.convolve()
  3. La integral de convolución con algoritmo

1. Un sistema con Scipy.signal.lti() y simulado con Scipy.signal.lsim()

Esta forma de describir el problema simplifica el desarrollar a partir de la descripción de la ecuación del sistema.

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación y condiciones de inicio en cero. La salida y(t) se calcula con una versión muestreada xi de la función de entrada x(t) a partir del tiempo muestreado en el intervalo [a,b].

En el algoritmo se incorpora la gráfica de la función de transferencia h(t) como referencia.

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)
import numpy as np
import scipy.signal as senal

# INGRESO
# tiempo [a,b]
a = 0 ; b = 5 ; muestras=101

# Señal de entrada
u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)

# Sistema
# coeficients Q  P de la ecuación diferencial
Q = [1., 3, 2]
P = [1., 0]
# condiciones de inicio [y'0,y0]
inicio = [0,0]

# PROCEDIMIENTO
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, inicio)
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.subplot(211)
plt.suptitle('Respuesta a estado cero')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')

plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(t_y, yi, color='magenta', label='y(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='y(t)')
# plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

2. La integral de convolución con Numpy.convolve()

El algoritmo requiere la expresión de las funciones x(t) y h(t) para determinar y(t). Realiza el cálculo de la convolución que para aproximar a la integral requiere multiplicar por dt.

Note que debe mantener el valor dt pequeño,  en el ejercicio se ha definido como 0.01 que corresponden a casi 1000 muestras.

Para usar la función convolve se requiere muestras en tiempos, simétricos alrededor de cero debido a la naturaleza del algoritmo. Sección que puede ser eliminada antes de mostrar la gráfica.  Para recordarlo se ha dejado la parte negativa del rango de tiempo.

# Lathi ejemplo 2.6 pdf129
# Señales contínuas en convolución
# Compare con el algoritmo en discretas
import numpy as np

# INGRESO
# tiempo [a,b] simétrico a 0
b = 5 ; a = -b; dt =0.01

u = lambda t: np.piecewise(t,t>=0,[1,0])

x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)


# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)

# Integral de Convolucion x[t]*h[t]
# corrección de magnitud por dt para en integral
yi = np.convolve(xi,hi,'same')*dt

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.figure(1)
plt.suptitle('Integral de Convolución x(t)*h(t)')

plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()

plt.show()

3. La integral de convolución con algoritmo

Esta sección desarrolla el algoritmo para la integral de convolución, tomando paso a paso las operaciones para obtener la gráfica del resultado. El algoritmo se presenta con propósitos didácticos.

Parece ser un poco más lento debido a que se realizan las operaciones de forma básica (interpretada) de Python, la función usa elementos compilados más eficientes.

Note que la única sección que cambia es la de convolución, por lo que la gráfica de salida es la misma que la sección anterior.:

# Lathi ejemplo 2.6 pdf129
# Señales contínuas en convolución
# Compare con el algoritmo en discretas
import numpy as np

# INGRESO
# Rango [a,b], simétrico a 0
b = 5 ; a = -b; dt =0.01

u = lambda t: np.piecewise(t,t>=0,[1,0])

x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)


# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)

# Integral de Convolución x[t]*h[t]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i in range(0,muestras):
    suma=0
    for k in range(0,muestras):
        suma = suma + x(ti[k])*h(ti[i]-ti[k])
    yi[i]= suma*dt

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.figure(1)
plt.suptitle('Integral de Convolución x(t)*h(t)')

plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()

plt.show()

3.2.1 Respuesta entrada cero – Python

1. Fórmulas simbólicas con Sympy

Para éste desarrollo se toma la función original , y se reordena de la forma:

(D^2 + 3D +2)y(t) = 0 D^2 y(t) + 3Dy(t) + 2y(t) = 0

Se define la variable t como símbolo, a y como función, el programa presenta la fórmula a la que se pueden aplicar las condiciones iniciales y una solución particular.

La ecuación diferencial: 
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
Solución general: 
           -t       -2*t
y(t) = C1*e   + C2*e    

Ecuaciones de condiciones iniciales
0 = C1 + C2
-5 = -C1 - 2*C2
{C1: -5, C2: 5}
Solución particular:
            -t      -2*t
y(t) = - 5*e   + 5*e    
>>>

Las intrucciones en Python son:

# Lathi 2.1.b pdf 116
# Respuesta a entrada cero con Sympy
# solucion para (D^2+ 3D + 2)y = 0
import sympy as sym

# INGRESO
# Condiciones iniciales
t0 = 0 ; y0 = 0 ; dy0 = -5

t = sym.Symbol('t')
y = sym.Function('y')
ecuacion = sym.Eq(sym.diff(y(t),t,t) + 3*sym.diff(y(t),t)+2*y(t),0)

# PROCEDIMIENTO
general = sym.dsolve(ecuacion, y(t))
general = general.expand()

# Condiciones iniciales
# rhs, righ hand side. Lado derecho de ecuacion
s0 = sym.Eq( y0, general.rhs.subs(t,t0))
s1 = sym.Eq(dy0, general.rhs.diff(t).subs(t,t0))
constante = sym.solve([s0,s1])

# Solucion particular
particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

# SALIDA
print('La ecuación diferencial: ')
sym.pprint(ecuacion)
print('Solución general: ')
sym.pprint(general)

print('\nEcuaciones de condiciones iniciales')
sym.pprint(s0)               
sym.pprint(s1)
print(constante)
print('Solución particular:')
sym.pprint(particular)

Para obtener la solución particular con las condiciones iniciales, se generan dos ecuaciones sutituyendo del lado izquierdo (lhs) los valores iniciales, y en el dado derecho (rhs) los valores de t=0 en la ecuación general y en la derivada del lado derecho.

Adicionalmente se puede realizar la gráfica, convirtiendo la ecuación de forma simbólica a forma numérica numpy con sym.lambdify(t,particular,’numpy’). Se establece el rango de observacion para evaluar la función y se envia a graficar. Se añaden las siguientes instrucciones:

# GRAFICA
import numpy as np
a = 0 ; b =5 ; muestras = 101
yt = sym.lambdify(t, particular.rhs, 'numpy')

ti = np.linspace(a,b,muestras)
yi = yt(ti)

import matplotlib.pyplot as plt
plt.plot(ti,yi, color='blue', label='y(t)')
plt.title('Respuesta entrada cero: y(t)='+str(particular.rhs))
plt.xlabel('t')
plt.grid()
plt.show()

2. Un sistema con Scipy.signal.lti() y simulado con Scipy.signal.lsim()

Esta forma de describir el problema simplifica el desarrollar a partir de la descripción de la ecuación del sistema.

La función Scipy.signal.lti() permite definir un sistema LTI, usando los coeficientes de la ecuación como numerador y denominador P(D)/Q(D).

La función scipy.signal.lsim(), usando:
– un intervalo [a,b] de tiempo ti muestreado,
– una señal de entrada x(t) como vector xi  y
– las condiciones de inicio  como vector [y'(0),y(0)]
determina la señal de salida en dos formas:
– la respuesta del sistema yi
– los componentes en forma de columnas yc=[yi(t),y0(t)]

usando los datos del ejercicio se grafica la respuesta a entrada cero:

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)
import numpy as np
import scipy.signal as senal


# INGRESO
# tiempo [a,b]
a = 0 ; b = 5 ; muestras=101

# Señal de entrada
u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: t*0 # 10*np.exp(-3*t)*u(t)

# Sistema
# coeficients Q  P de la ecuación diferencial
Q = [1., 3, 2]
P = [1., 0]
# condiciones de inicio [y'0,y0]
inicio = [-5,0]

# PROCEDIMIENTO
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, inicio)
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.subplot(211)
plt.suptitle('Respuesta a Entrada cero')
plt.plot(ti, xi, color='blue', label='x(t)')
#plt.plot(t_h, hi, color='red', label='h(t)')

plt.legend()
plt.grid()

plt.subplot(212)
# plt.plot(t_y, yi, color='magenta', label='y(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='y(t)')
plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

3. Algoritmo de Runge-Kutta d2y/dx2 de Análisis numérico

Se reordena el sistema de ecuaciones para plantear una solución con el algoritmo Runge-Kutta de 2do orden para ecuaciones diferenciales con segundas derivadas

(D^2 + 3D +2)y(t) = 0 D^2 y(t) + 3Dy(t) + 2y(t) = 0

Se ubica el término de mayor grado a la izquierda de la ecuación:

D^2 y(t) = - 3Dy(t) - 2y(t)

derivadas sutituyendo:

Dy = y’ = z = f(x)
D2y = y» = z’= -3z – 2y = g(x)

con los valores iniciales de t0=0, y0=0, z0=-5 evaluado en el rango de tiempo [0,tn=5] y como tamaño de paso h=dt=0.1

El algoritmo permite obtener la gráfica y la tabla de datos.

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = -5

h = 0.1
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(ti,yi, label='y(t)')

plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Más detalles del algoritmo en :

8.2.2 Runge-Kutta d2y/dx2

3.2.2 Respuesta entrada cero – Diagrama Bloques

Desarrollo con diagrama de bloques y un simulador

Un diferenciador D representa un sistema inestable pues para una entrada acotada como el escalón μ(t) resulta en una salida no acotada como el impulso unitario δ(t). Un sistema con diferenciadores expuesto a ruido amplifica el ruido.

Para los diagramas de bloques y simuladores se prefiere usar integradores (1/D). Con problema planteado ésta sección se cambia la forma de la ecuación:

D^2 y(t) + 3Dy(t) + 2y(t) = 0

Se ubica el término de mayor grado a la izquierda de la ecuación:

D^2 y(t) = - 3Dy(t) - 2y(t)

se despeja la respuesta y(t):

\frac{D^2 y(t)}{D^2} = \frac{-3Dy(t) - 2y(t)}{D^2} y(t) = -3 \frac{1}{D} y(t) - 2 \frac{1}{D^2} y(t)

La forma de la ecuación permite crear el diagrama del sistema, usando 1/D como un integrador y 1/D2 representa una segunda integración en serie.

Los coeficientes de los términos se ubican como ganancias retroalimentadas. Los valores iniciales se configuran en los integradores observando cual salida se usará en el sumador.

Usando el Simulador XCOS-Scilab

Para XCOS, se añade un visor de señal de salida (Sinks/Scope) y un reloj de muestreo. Para y'(0)=5 siguiendo el diagrama se ubica y'(t)=dy/dt como la salida del primer integrador . El valor se cambia dando doble click al elemento marcado en rojo y escribiendo el valor de -5.

La gráfica se obtiene usando el osciloscopio (scope)

El tiempo de observación se establece con el valor de «Simulatión/Final Integration time» a 5 segundos.


De forma semejante se puede realizar el ejercicio en otros simuladores, como por ejemplo Simulink-Matlab.


4. Desarrollo como circuito en simulador

La simulación en forma de circuito usando las herramientas de Simulink/Matlab y SimPowerSystems, para cada componente (Elements) se crea con «Series RLC Branch». Para extraer las lecturas de las señal se usa un medidor (Measurements) que se ubica en serie con los demás componentes y se visializa con un osciloscopio (Simulink/Commonly Used Blocks/Scope).

Las conexiones son similares a las que se usarían en el laboratorio. Se puede observar que las respuestas serán las mismas.

Las condiciones iniciales se pueden establecer en los componentes que pueden tener un valor inicial como el capacitor en 5 voltios. El signo de la corriente de salida ya está incluido en la inversión de polaridad respecto a la polaridad de la fuente.

3.1 Sistema – Modelo entrada-salida

Referencia: Lathi 2.1 pdf114, Schaum/Hsu 2.5.B p60

La descripción de un sistema en términos de mediciones en los extremos se denomina Modelo de entrada-salida

.

Para el caso de un circuito eléctrico RLC, la creación del modelo inicia con la descripción de la ecuación que relaciona el voltaje x(t) de entrada y la corriente de salida y(t).

La respuesta de un sistema lineal puede ser expresado como la suma de dos componentes: respuesta a entrada cero y respuesta a estado cero.

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

Ejemplo

Referencia: Lathi 1.8-1 pdf/p.80. Oppenheim problema 2.61c pdf/p.191

Como ejemplo, usaremos la corriente de lazo del circuito mostrado.

FIEC05058_RLCusando la ley de voltajes en lazo se tiene que:

vL(t) + vR(t) +vC(t) = x(t)

que con las leyes de corriente para cada elemento (inductor, resistor y capacitor) se traducen en:

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

Para tener todo en función de un solo operador, se derivan ambos lados de la ecuación:

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

que es la relación de «entrada-salida» del sistema con la entrada x(t) y la salida y(t) que permitirá analizar el sistema que representa el circuito.


Notación D y 1/D

Por conveniencia, para usar una notación más compacta de la ecuación diferencial, el operador dy/dt se cambia por la notación D.

\frac{dy}{dt} = Dy(t) \frac{d^2 y}{dt^2} = D^{2}y(t)

Que convierte la ecuación de entrada-salida a la expresión:

(D^2 + 3D +2)y(t) = Dx(t)

que es identica a la expresión de entrada y salida que describe al circuito.

El operador diferencial es posible usarlo para la notación con integrales,

\int y(t)dt = \frac{1}{D}y(t)

por lo que la expresión del circuito también se puede escribir como

\Big( D + 3D + \frac{2}{D} \Big) y(t) = x(t)

al multiplicar ambos lados por el operador D se convierte nuevamente en una expresión sin denominadores D.

Recuerde: la expresión con operadores D, NO ES una ecuación algebraica, pues la expresión de operadores D aplican solo a y(t).

En adelante, el sistema descrito por ecuaciones diferenciales usa el operador D=\frac{d}{dt}, por ejemplo:

a_2 \frac{d^2}{dt^2}y(t) + a_1 \frac{d}{dt}y(t) + a_0 y(t) = b_1\frac{d}{dt}x(t) + b_0x(t) a_2 D^ 2y(t) + a_1 Dy(t) + a_0y(t) = b_1Dx(t) + b_0 x(t) (a_2 D^ 2 + a_1 D + a_0)y(t) = (b_1D + b_0 )x(t) Q(D) y(t) = P(D) x(t)

3.2 Respuesta entrada cero-LTIC

Referencia: Lathi 2.1 pdf114, Schaum/Hsu 2.5 p60

Siguiendo el modelo entrada-salida del sistema, la respuesta se obtiene como la superposición de sus componentes:

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

Para entrada cero se usa x(t)=0, es decir al sistema no se le aplica una señal de entrada, se pone a tierra la entrada y se observa la salida. El sistema, al no tener señal entrada x(t) simplifica el problema a:

Q(D) y(t) = P(D) x(t) Q(D) y(t) = 0

Ejemplo

Referencia: Lathi Ejemplo 2.1.a pdf116, Ejemplo 2.2 pdf 119

Encuentre la respuesta a entrada cero para el sistema LTIC del circuito, descrito por la ecuación:

(D^2 + 3D +2)y(t) = Dx(t)

con las condiciones iniciales de y0(t) =0 , y’0(t) =-5
La entrada cero, x(t)=0 convierte el circuito en:

Se plantean varias formas de desarrollo, para explorar algunas opciones usando los conceptos descritos en la parte teórica:

  1. Desarrollo analítico
  2. Desarrollo usando Python:
    1. con fórmulas simbólicas con Sympy
    2. con Runge-Kutta d2y/dx2 de análisis numérico
  3. Desarrollo usando simulador con:
    1. con diagrama de bloques
    2. como circuito en simulador

1. Desarrollo analítico

Siendo el sistema descrito por:

Q(D) y(t) = P(D) x(t)

Para entrada cero se usa x(t)=0, se quita la fuente y se cierra el circuito, observando solo lo que hace el sistema sin señal de entrada. De ésta forma se simplifica la ecuación al eliminar P(D):

(D^2 + 3D +2)y(t) = 0

El polinomio Q(D) es el polinomio característico del sistema, entonces la ecuación característica del sistema es:

λ2 + 3 λ + 2 = (λ+1)(λ+2)
Raíces características Modos característicos
λ1 = -1 e-t
λ2 = -2 e-2t

con los resultados característicos se forma la respuesta a entrada cero como y0(t).

y_0 (t) = c_1 e^{-t} + c_2 e^{-2t}

Para determinar las constantes se usan las condiciones iniciales. No hay señal en la salida al tiempo cero y(0)=0, y para observar que sucedería si el circuito tuviese energía interna, por ejemplo un capacitor cargado se tiene que y'(0)=-5:

y_0 (t) = c_1 e^{-t} + c_2 e^{-2t} y'_0 (t) = -c_1 e^{-t} -2c_2 e^{-2t}
evaluando en condiciones iniciales:
 0 =  c1 +  c2
-5 = -c1 - 2c2

se resuelve obeniendo:
c1 = -5 
c2 = 5

que al sustituir en la ecuacion inicial, se tiene la respuesta a entrada cero:

y_0(t) = -5e^{-t} +5e^{-2t}

La gráfica muestra el sentido de la corriente, usando la carga residual del capacitor dentro del circuito, a pesar de no tener señal de entrada a partir del tiempo 0 hasta 5:

Para determinar los valores de las constantes, se puede usar algunas instrucciones sencillas en Pyton y Numpy, para las raices del polinomio se usan solo los coeficioentes de grado mayor a menor.

>>> import nupmpy as np

>>> np.roots([1,3,2])
array([-2., -1.])

Para los coeficientes se plantean las ecuaciones de la forma matricial Ax=B:

>>> A = [[ 1, 1],
	 [-1,-2]]
>>> B =  [ 0,-5]
>>> np.linalg.solve(A,B)
array([-5.,  5.])
>>>

La parte de las gráficas se desarrollan en las secciones dedicadas usar Python para obtener las respuestas.

3.4 Respuesta a estado cero – LTIC

Referencia: Lathi 2.4 pdf124, Oppenheim 2.1 p97/pdf125 , Schaum/Hsu 2.5.B p.60

El estado cero del sistema, «Zero-State», supone no hay energía almacenada, que los capacitores están descargados, que recien sale el equipo de la caja.  Para éste caso, la respuesta del sistema se conoce como respuesta a estado cero, «Zero-State response».

Para los problemas presentados se asume que el sistema es lineal , causal e invariante en el tiempo. En la práctica, muchos de los sistemas son causales, pues su respuesta no inicia antes de aplicar una entrada, es decir, todas las entradas a evaluar empiezan en t=0.

la respuesta del sistema y(t) para un LTIC se determina como:

y(t) = x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) d\tau

Es importante observar que el integral de convolución se realiza con respecto a τ en lugar de t.

Si la entrada x(t) y el sistema h(t) son causales, la respuesta también será causal.

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) d\tau , & t\ge 0\\ 0, & t<0 \end{cases}

El límite inferior del integral se usa como 0, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.


Ejemplo 01

Referencia:  Lathi ejemplo 2.6 pdf129

Encuentre la corriente y(t) del circuito RLC, cuando todas las condiciones iniciales son cero y en la entrada se tiene la señal x(t) descrita por:

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

Además el sistema tiene respuesta a impulso:

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

El ejemplo es la continuación del presentado para respuesta a entrada cero, que tiene la ecuación:

(D^2 + 3D +2)y(t) = Dx(t)

Desarrollo Analítico

La respuesta se obtiene aplicando convolución:

y(t) = x(t) \circledast h(t) = [ 10 e^{-3t} \mu (t)] \circledast [(2e^{-2t} - e^{-t}) \mu (t)]

usando la propiedad distributiva de la convolución:

y(t) = [10e^{-3t} \mu (t) \circledast 2e^{-2t} \mu (t)] - [10e^{-3t} \mu (t) \circledast e^{-t} \mu (t)] = 20[e^{-3t}\mu (t) \circledast e^{-2t} \mu (t)] - 10[e^{-3t} \mu(t) \circledast e^{-t} \mu (t)]

Para éste ejercicio se usa la tabla de convoluciones, asi se enfoca en la forma de la señal resultante, el siguiente ejemplo se desarrolla el integral de convolución.

y(t) = 20\frac{e^{-3t} - e^{-2t}}{-3-(-2)}\mu (t) - 10\frac{e^{-3t} - e^{-t}}{-3-(-1)}\mu (t) = -20[e^{-3t} - e^{-2t}]\mu (t) + 5[e^{-3t} - e^{-t}]\mu (t) = [-5e^{-t} + 20e^{-2t} - 15e^{-3t}]\mu (t)

De las gráficas se observa que la entrada es semejante a conectar en la entrada un capacitor con carga, que la pierde en el tiempo.

En la salida se observa el efecto, la parte inicial corresponde a la corriente en el circuito mientras el capacitor de la entrada entrega energía al sistema. Note que en el sistema o circuito se debe ir cargando el capacitor del sistema. Luego, un poco más del segundo 1, la corriente invierte el sentido volviéndose negativa por la carga almacenada en el capacitor del sistema.

La explición breve realizada debería ser comprobada en los experimentos de laboratorio, preferiblemente a escala menor con componentes tipo electrónico.
Proponer como tarea.

Ejemplo 02

Referencia: Lathi Ejemplo 2.5 pdf127

Para un sistema LTIC, si la respuesta al impulso es

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

determine la respuesta y(t) para la entrada

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

Desarrollo Analítico

Para éste ejercicio se desarrollará la integral de convolución. La entrada y respuesta al impulso se convierte a:

x(\tau) = e^{-\tau} u(\tau) h(t-\tau) = e^{-2(t-\tau)} u(t-\tau)

recuerde que la integración es respecto a τ enel intervalo 0≤τ≤t.

y(t) = \begin{cases} \int_{0}^{t} e^{-\tau}u(\tau) e^{-2(t-\tau)}u(t-\tau) d\tau , & t\ge 0 \\0, & t \lt 0 \end{cases}

Los valores de u(τ) =1 debido a se convierte a 0 para τ<0 y en el caso de u(t-τ)=1 se conviertea 0 cuando τ≥t.

y(t) = \int_{0}^{t}e^{-\tau} e^{-2(t-\tau)} d\tau = \int_{0}^{t} e^{-\tau} e^{-2t} e^{2\tau} d\tau = e^{-2t} \int_{0}^{t} e^{\tau} d\tau = e^{-2t} e^{\tau} \Big|_0^t

Evaluado en el rango de integración:

= e^{-2t} (e^{t} - 1) = e^{-t} - e^{-2t}

para t≥0, y además como y(t) = 0 para t<0

y(t) = (e^{-t} - e^{-2t})u(t)

3.5 Respuesta a impulsos – discreto

Referencia: Lathi 3.8 pdf193, Oppenheim 2.1.2 p77/pdf105, Schaum 2.6.C p62,

Una señal discreta de entrada x[n] se representa como una superposición de versiones escaladas de un conjunto de impulsos unitarios desplazados δ[n-k], cada uno con valor diferente de cero en un solo punto en el tiempo, especificado por el valor de k.

La respuesta de un sistema lineal y[n] a x[n] es la superposición de las respuestas escaladas del sistema a cada uno de estos impulsos desplazados.

Si usamos hk[n] como la respuesta del sistema lineal al impulso unitario desplazado por δ[n-k]. La respuesta y[n] del sistema lineal a la entrada x[n] en la ecuacion será la combinación lineal ponderada de las respuesta básicas.

y[n] = \sum_{k=-\infty}^{+\infty} x[k]h_{k}[n]

Si se conoce la respuesta de un sistema lineal al conjunto de impulsos unitarios desplazados, podemos construir una respuesta a una entrada arbitraria.

Si el sistema lineal también es invariante en el tiempo, entonces estas respuestas a impulsos unitarios desplazados en el tiempo son todas las versiones desplazadas en el tiempo unas de otras.

h[n] es la salida del sistem LTI cuando δ[n] es la entrada. Entonces para un sistema LTI la ecuación se vuelve.

y[n] = \sum_{k=-\infty}^{+\infty} x[k]h[n-k]

El resultado se conoce como la «suma de convolución» o «suma de superposición» y la operación miembro derecho de la ecuación se llama convolución de las secuencias x[n] y h[n] que se representa de manera simbólica como:

y[n] = x[n]*h[n]

Ejemplo

Referencia: Openheim Ejemplo 2.4 p85/pdf113

Considere una entrada x[n] y una respuesta al impulso unitario h[n] dada por:

x[n] = (u[n]-u[n-5])
h[n] = αn (u[n]-u[n-7])

con α>1.  Para el ejercicio, α=1.5.

Nota: si α es un entero, por ejemplo 2, usar α=2.0, para que la operación potencia se realice con números reales, como entero se puede saturar y dar error.


Al aplicar la suma de convolución obtiene:

Para aplicar el algoritmo se requiere definir u[n], por ser parte de las funciones x[n] y h[n]. Dado que la operación requiere valores fuera del rango muestreado para n, la sección suma convolución utiliza las funciones en lugar de los vectores xi, hi.

La función está definida en un intervalo simétrico, por lo que el rango de trabajo [a,b] se mantiene de la forma [-b,b] en las intrucciones.

# Ejemplo Oppenheim Ejemplo 2.3 p83/pdf111
import numpy as np

# INGRESO
# Rango [a,b], simétrico a 0
b = 15 ; a = -b

alfa = 1.5
u = lambda n: np.piecewise(n,n>=0,[1,0])

x = lambda n: u(n)-u(n-5)
h = lambda n: (alfa**n)*(u(n)-u(n-7))


# PROCEDIMIENTO
ni = np.arange(a,b+1,1)
xi = x(ni)
hi = h(ni)

# Suma de Convolucion x[n]*h[n]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i in range(0,muestras):
    suma=0
    for k in range(0,muestras):
        suma = suma + x(ni[k])*h(ni[i]-ni[k])
    yi[i]= suma

# yi = np.convolve(xi,hi,'same')

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.figure(1)
plt.suptitle('Suma de Convolución x[n]*h[n]')

plt.subplot(311)
plt.stem(ni,xi,
         linefmt='b--',
         markerfmt='bo',
         basefmt='k-')
plt.ylabel('x[n]')

plt.subplot(312)
plt.stem(ni,hi,
         linefmt='b--',
         markerfmt='ro',
         basefmt='k-')
plt.ylabel('h[n]')

plt.subplot(313)
plt.stem(ni,yi,
         linefmt='g-.',
         markerfmt='mo',
         basefmt='k-')
plt.ylabel('x[n]*h[n]')
plt.xlabel('n')

plt.show()

La suma convolución se encuentra también disponible con Numpy en np.convolve(), la  sección de suma convoluacion se puede reemplazar y obtener los mismos resultados. Considere que para éste caso se usan los vectores xi y hi.

yi = np.convolve(xi,hi,'same')

el algoritmo se puede aplicar a otros ejercicio para comprobar los resultados.

Continúe con el ejercicio 2.5 del libro.