Uso de ipython para el modelado de un sistema a partir de los datos obtenidos en un ensayo. Ensayo en lazo abierto
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 [4]:
#Abrimos el fichero csv con los datos de la muestra para filtrarlos
err_u = 2
err_d = 1
datos_sin_filtrar = pd.read_csv('datos.csv')
datos = datos_sin_filtrar[(datos_sin_filtrar['Diametro X'] >= err_d) & (datos_sin_filtrar['Diametro Y'] >= err_d) & (datos_sin_filtrar['Diametro X'] <= err_u) & (datos_sin_filtrar['Diametro Y'] <= err_u)]
datos.describe()
Out[4]:
In [5]:
#Almacenamos en una lista las columnas del fichero con las que vamos a trabajar
#columns = ['temperatura', 'entrada']
columns = ['Diametro X', 'RPM TRAC']
Representamos los datos mostrados en función del tiempo. De esta manera, vemos la respuesta física que tiene nuestro sistema.
In [6]:
#Mostramos en varias gráficas la información obtenida tras el ensayo
th_u = 1.95
th_d = 1.55
datos[columns].plot(secondary_y=['RPM TRAC'],figsize=(10,5),title='Modelo matemático del sistema').hlines([th_d ,th_u],0,2000,colors='r')
#datos_filtrados['RPM TRAC'].plot(secondary_y=True,style='g',figsize=(20,20)).set_ylabel=('RPM')
Out[6]:
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 [7]:
# Buscamos el polinomio de orden 4 que determina la distribución de los datos
reg = np.polyfit(datos['time'],datos['Diametro X'],4)
# Calculamos los valores de y con la regresión
ry = np.polyval(reg,datos['time'])
print (reg)
In [8]:
d = {'Diametro X' : datos['Diametro X'],
'Ajuste': ry, 'RPM TRAC' : datos['RPM TRAC']}
df = pd.DataFrame(d)
df.plot(subplots=True,figsize=(20,20))
Out[8]:
In [9]:
plt.figure(figsize=(10,10))
plt.plot(datos['time'],datos['Diametro X'], 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[9]:
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 [10]:
##Respuesta en frecuencia del sistema
In [11]:
num = [25.9459 ,0.00015733 ,0.00000000818174]
den = [1,0,0]
In [12]:
tf = signal.lti(num,den)
w, mag, phase = signal.bode(tf)
In [13]:
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[13]:
In [14]:
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[14]:
In [15]:
t, y = signal.step2(tf) # Respuesta a escalón unitario
plt.plot(t, 2250 * y) # Equivalente a una entrada de altura 2250
Out[15]:
In [ ]: