Modelado de un sistema con ipython

Uso de ipython para el modelado de un sistema a partir de los datos obtenidos en un ensayo. Ensayo en lazo abierto


In [22]:
#Importamos las librerías utilizadas
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import signal

In [23]:
#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__))


Numpy v1.9.2
Pandas v0.16.2
Seaborn v0.6.0

In [24]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [62]:
#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[62]:
time Tmp Husillo Tmp Nozzle Diametro X Diametro Y MARCHA PARO RPM EXTR RPM TRAC
count 1299.000000 1299.000000 1299.000000 1299.000000 1299.000000 1299 1299 1299.000000 1299.000000
mean 1110.278676 64.473826 138.144419 1.570151 1.503133 1 1 1.983834 3.575253
std 552.645419 0.354294 0.806497 0.188291 0.186580 0 0 0.126163 0.789455
min 235.000000 63.500000 135.600000 1.000410 1.000236 True True 1.000000 3.070136
25% 576.500000 64.300000 137.600000 1.470675 1.402492 1 1 2.000000 3.219072
50% 1134.000000 64.400000 138.100000 1.585374 1.494436 1 1 2.000000 3.219072
75% 1622.500000 64.800000 138.800000 1.694338 1.632352 1 1 2.000000 3.427583
max 2003.000000 65.100000 140.300000 1.998289 1.988636 True True 2.000000 5.899920

In [40]:
#Almacenamos en una lista las columnas del fichero con las que vamos a trabajar
#columns = ['temperatura', 'entrada']
columns = ['Diametro X', 'RPM TRAC']

Representación

Representamos los datos mostrados en función del tiempo. De esta manera, vemos la respuesta física que tiene nuestro sistema.


In [64]:
#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[64]:
<matplotlib.collections.LineCollection at 0x149b6510>

Cálculo del polinomio característico

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 [71]:
# 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)


[ -2.54072222e-13   1.07702622e-09  -1.48292474e-06   7.58672664e-04
   1.43592440e+00]

In [94]:
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[94]:
array([<matplotlib.axes._subplots.AxesSubplot object at 0x17536510>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x1756CFD0>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x1759C850>], dtype=object)

In [86]:
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[86]:
<matplotlib.text.Text at 0x16988e50>

El polinomio caracteristico de nuestro sistema es: $P_x= 25.9459 -1.5733·10^{-4}·X - 8.18174·10^{-9}·X^2$

Transformada de laplace

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 [ ]:
num = [25.9459 ,0.00015733 ,0.00000000818174]
den = [1,0,0]

In [ ]:
tf = signal.lti(num,den)
w, mag, phase = signal.bode(tf)

In [ ]:
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

In [ ]:
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')

In [ ]:
t, y = signal.step2(tf) # Respuesta a escalón unitario
plt.plot(t, 2250 * y) # Equivalente a una entrada de altura 2250

In [ ]: