Modelado de un sistema con ipython

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__))


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

In [3]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

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']

Representación

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]:
<matplotlib.axes._subplots.AxesSubplot at 0x1336f1b0>

Cálculo del polinomio

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)


[  8.18174358e-09  -1.57332539e-04   2.59458703e+01]

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]:
<matplotlib.text.Text at 0x13eeb070>

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 [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]:
[<matplotlib.lines.Line2D at 0x153fc350>]

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]:
[<matplotlib.lines.Line2D at 0x142e71b0>]

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]:
[<matplotlib.lines.Line2D at 0x1421fdb0>]

In [ ]: