<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Simulacion PID</span> por <a xmlns:cc="http://creativecommons.org/ns#" href="http://nbviewer.ipython.org/urls/raw.githubusercontent.com/elsuizo/Redes_neuronales_Fuzzy/master/guia1.ipynb?create=1" property="cc:attributionName" rel="cc:attributionURL">Martin Noblía Marcos Panizzo Damián Presti</a> se distribuye bajo una Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.
Vamos a implementar la simulación numérica de un controlador PID siguiendo el siguiente articulo:
Discretization of simulator, filter, and PID controller. Finn Haugen (http://techteach.no/publications/articles/discretization/discretization.pdf)
In [76]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Un sistema de segundo orden en su forma estándar se define:
$$\frac{Y(s)}{R(s)}=\frac{\omega_{n}^{2}}{s^{2}+2\,\zeta\omega_{n}s + \omega_{n}^{2}}$$Luego el comportamiento dinámico del sistema de segundo orden se describe en términos de dos parámetros $\zeta$ y $\omega_{n}$.
Si $0< \zeta < 1$, los polos a lazo cerrado son complejos conjugados y se encuentran en el semiplano izquierdo del plano $s$ (BIBO-estable)
y los siguientes casos:
subamortiguado: $0< \zeta < 1$
criticamente amortiguado: $\zeta = 1$
sobreamortiguado: $\zeta > 1$
Sabemos que la salida de un sistema de segundo orden en el espacio temporal es:
$$y(t) =1 - e^{-\zeta \omega_{n}t}(cos(\omega_{d} \, t)+\frac{\zeta}{\sqrt{1-\zeta^{2}}}sin(\omega_{d} \, t)) $$La señal de error se define como la diferencia entre la entrada y la salida:
$$e(t) = r(t) - y(t) = e^{-\zeta \omega_{n}t}(cos(\omega_{d} \, t)+\frac{\zeta}{\sqrt{1-\zeta^{2}}}sin(\omega_{d} \, t))$$A continuación definimos estas señales como funciones que dependen de los parámetros caracteristicos que vimos.
In [77]:
def y(t, zeta, omega_d):
return 1 - np.exp(-zeta * omega_d * t) * (np.cos(omega_d * t) + (zeta / np.sqrt(1 - zeta**2)) * np.sin(omega_d * t))
In [78]:
def error(t, zeta, omega_d):
return np.exp(-zeta * omega_d * t) * (np.cos(omega_d * t) + (zeta / np.sqrt(1 - zeta**2)) * np.sin(omega_d * t))
Siguiendo el articulo sabemos que finalmente la discretización del controlador PID queda:
$$u(t_{k})=u(t_{k-1})+[u_{0}(t_{k})-u_{0}(t_{k-1})]+K_{p}[e(t_{k})-e(t_{k-1})]+\frac{K_{p}T_{s}}{T_{i}}e(t_{k})+\frac{K_{p}T_{d}}{T_{s}}[e(t_{k})-2e(t_{k-1})+e(t_{k-2})] $$Luego el controlador total o incremental se calcula:
$$\Delta u(t_{k})=[u_{0}(t_{k})-u_{0}(t_{k-1})]+K_{p}[e(t_{k})-e(t_{k-1})]+ \frac{K_{p} T_{s}}{T_{i}}e(t_{k})+\frac{K_{p}T_{d}}{T_{s}}[e(t_{k})-2e(t_{k-1})+e(t_{k-2})]$$$$u(t_{k})=u(t_{k-1})+\Delta u(t_{k})$$A continuación la implementación:
In [79]:
N = 1000 # numero de muestras
T_s = .1 # Periodo de muestreo
k = np.arange(N) # vector discreto
u = np.zeros(N) # vector de la ley de control
K_p = 1 # Constante proporcional
T_i = .73 # Constante de tiempo de integracion
T_d = .1 # Constante de actualizacion de la derivada
t_k = T_s * k # Tiempo discretizado
In [80]:
D_u = np.zeros_like(u)
# Parametros del sistema de segundo orden
zeta = .3
omega = 0.9
# Simulacion
for k in xrange(0,len(t_k)):
D_u[k] = (K_p * (error(t_k[k], zeta, omega) - error(t_k[k-1], zeta, omega)) + ((K_p * T_s) / T_i) * error(t_k[k],zeta, omega) +
((K_p * T_d) / T_s) * (error(t_k[k], zeta, omega) - 2 * error(t_k[k-1],zeta,omega) + error(t_k[k-2],zeta,omega)) )
u[k] = u[k-1] + D_u[k]
In [81]:
# Graficos
plt.rcParams['figure.figsize'] = 8,6 #parámetros de tamaño
plt.plot(t_k, u)
plt.plot(t_k, y(t_k, zeta, omega))
plt.plot(t_k, error(t_k, zeta, omega))
plt.xlabel('$t$', fontsize=20)
plt.title('Simulacion de un Controlador PID', fontsize=20)
plt.legend(['$u(t)$','$y(t)$', '$e(t)$'],prop={'size':20})
plt.grid()
In [89]:
# Animacion de la señal y(t)
# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation
fig = plt.figure()
ax = plt.axes(xlim=(0, 100), ylim=(-2, 2))
line, = ax.plot([], [], lw=2)
def init():
line.set_data([], [])
return line,
# variamos el parametro
def animate(i):
line.set_data(t_k, y(t_k, zeta/(i+1), omega ))
return line,
animation.FuncAnimation(fig, animate, init_func=init,
frames=100, interval=20, blit=True)
Out[89]:
In [82]: