Uso de ipython para el modelado de un sistema a partir de los datos obtenidos en un ensayo.
In [1]:
#Importamos las librerías utilizadas
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import signal
In [2]:
#Mostramos las versiones usadas de cada librerías
print ("Numpy v{}".format(np.__version__))
print ("Pandas v{}".format(pd.__version__))
print ("Seaborn v{}".format(sns.__version__))
In [3]:
%pylab inline
In [6]:
#Abrimos el fichero csv con los datos de la muestra
datos = pd.read_csv('datos.csv')
In [7]:
#Almacenamos en una lista las columnas del fichero con las que vamos a trabajar
#columns = ['temperatura', 'entrada']
columns = ['temperatura', 'entrada']
Representamos los datos mostrados en función del tiempo. De esta manera, vemos la respuesta física que tiene nuestro sistema.
In [24]:
#Mostramos en varias gráficas la información obtenida tras el ensayo
datos[columns].plot(secondary_y=['entrada'],figsize=(10,5), ylim=(20,60),title='Modelo matemático del sistema')
#datos_filtrados['RPM TRAC'].plot(secondary_y=True,style='g',figsize=(20,20)).set_ylabel=('RPM')
Out[24]:
Hacemos una regresión con un polinomio de orden 4 para calcular cual es la mejor ecuación que se ajusta a la tendencia de nuestros datos.
In [34]:
# Buscamos el polinomio de orden 4 que determina la distribución de los datos
reg = np.polyfit(datos['time'],datos['temperatura'],2)
# Calculamos los valores de y con la regresión
ry = np.polyval(reg,datos['time'])
print (reg)
In [28]:
plt.plot(datos['time'],datos['temperatura'],'b^', label=('f(x)'))
plt.plot(datos['time'],ry,'ro', label=('regression'))
plt.legend(loc=0)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
Out[28]:
El polinomio caracteristico de nuestro sistema es: $P_x= 25.9459 -1.5733·10^{-4}·X - 8.18174·10^{-9}·X^2$
Si calculamos la transformada de laplace del sistema, obtenemos el siguiente resultado: \newline $G_s = \frac{25.95·S^2 - 0.00015733·S + 1.63635·10^{-8}}{S^3}$
In [ ]:
##Respuesta en frecuencia del sistema
In [37]:
num = [25.9459 ,0.00015733 ,0.00000000818174]
den = [1,0,0]
In [44]:
tf = signal.lti(num,den)
w, mag, phase = signal.bode(tf)
In [45]:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 6))
ax1.semilogx(w, mag) # Eje x logarítmico
ax2.semilogx(w, phase) # Eje x logarítmico
Out[45]:
In [49]:
w, H = signal.freqresp(tf)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 10))
ax1.plot(H.real, H.imag)
ax1.plot(H.real, -H.imag)
ax2.plot(tf.zeros.real, tf.zeros.imag, 'o')
ax2.plot(tf.poles.real, tf.poles.imag, 'x')
Out[49]:
In [50]:
t, y = signal.step2(tf) # Respuesta a escalón unitario
plt.plot(t, 2250 * y) # Equivalente a una entrada de altura 2250
Out[50]:
In [ ]: