En esta práctica nos enfocaremos a entender las metodologías para hacer y sintonizar controladores, para lo cual necesitaremos empezar con sistemas simples, nuestro primer ejemplo será la ecuación diferencial:
$$ \dot{x}(t) + x(t) = u(t) $$Lo primero que hacemos es aplicar una transformada de Laplace para obtener la siguiente función de transferencia:
$$ G(s) = \frac{X(s)}{U(s)} = \frac{1}{s+1} $$Una vez que tenemos esta función de transferencia, podemos ingresarla a una función de transferencia por medio de la librería control, de donde obtenemos la función tf:
In [ ]:
from control import tf
Para utilizar esta función, le damos un arreglo con los coeficientes del numerador y denominador de la función de transferencia, de tal manera que:
In [ ]:
G1 = tf([1],[1,1])
G1
Ya que tenemos definida la función de transferencia, podemos simular su comportamiento por medio de la función step_response:
In [ ]:
from control import step_response
In [ ]:
t, y = step_response(G1)
Y una vez simulado, graficar los datos:
In [ ]:
%matplotlib inline
In [ ]:
from matplotlib.pyplot import figure, close, style
style.use("ggplot")
In [ ]:
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y)
Sin embargo un sistema tan sencillo no tiene el comportamiento que necesitamos para este analisis; podemos utilizar el sistema masa resorte amortiguador para el siguiente ejemplo, recordemos que su función de transferencia es:
$$ G(s) = \frac{1}{ms^2 + cs + k} $$por lo que definimos su función de tranferencia, simulamos y graficamos:
In [ ]:
m = 1000
c = 1500
k = 15000
G2 = tf([1], [m, c, k])
In [ ]:
t, y = step_response(G2)
In [ ]:
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y)
Basado en la función de transferencia de un sistema masa resorte amortiguador:
In [ ]:
# ESCRIBE TU CODIGO AQUI
raise NotImplementedError
t, y = step_response(G3)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y);
In [ ]:
from numpy import sqrt
num = G3.num[0][0]/G3.den[0][0][0]
denom = (G3.den[0][0])/G3.den[0][0][0]
assert len(denom) is 3, "¿Estas usando la función del sistema masa resorte amortiguador?"
ωn = sqrt(denom[2])
ζ = denom[1]/(2*ωn)
assert 0 < ζ < 1, "Este sistema no es subamortiguado"
In [ ]:
# ESCRIBE TU CODIGO AQUI
raise NotImplementedError
t, y = step_response(G4)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y);
In [ ]:
from numpy import sqrt
num = G4.num[0][0]/G4.den[0][0][0]
denom = (G4.den[0][0])/G4.den[0][0][0]
assert len(denom) is 3, "¿Estas usando la función del sistema masa resorte amortiguador?"
ωn = sqrt(denom[2])
ζ = denom[1]/(2*ωn)
assert ζ > 1, "Este sistema no es sobreamortiguado"
In [ ]:
# ESCRIBE TU CODIGO AQUI
raise NotImplementedError
t, y = step_response(G5)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y);
In [ ]:
from numpy.testing import assert_array_equal
from numpy import array
polos = G5.pole()
polos_real = [p.real for p in polos]
polos_crit = array([p == 0 for p in polos_real])
polos_ines = array([p > 0 for p in polos_real])
assert polos_crit.any() and not polos_ines.any(), "Este sistema no es criticamente estable"
In [ ]:
# ESCRIBE TU CODIGO AQUI
raise NotImplementedError
t, y = step_response(G5)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y);
In [ ]:
from numpy.testing import assert_array_equal
from numpy import array
polos = G6.pole()
polos_real = [p.real for p in polos]
polos_crit = array([p == 0 for p in polos_real])
polos_ines = array([p > 0 for p in polos_real])
assert polos_ines.any(), "Este sistema no es inestable"
Podemos incluso graficar el comportamiento del sistema con diferentes parametros al mismo tiempo; para esto definiremos varios sistemas con una list comprehension:
In [ ]:
m = 1000
k = 15000
cs = [1500, 500, 0, -100, -500]
F = 1000
Gs = [tf([F], [m,c,k]) for c in cs]
Para asegurar que todas las simulaciones se harán con los mismos tiempos iniciales y finales, definiremos un arreglo desde $0$ hasta $10$ con $1000$ puntos para simular:
In [ ]:
from numpy import linspace
ts = linspace(0, 10, 1000)
Entonces, podremos simular el comportamiento de cada uno de estos sistemas con otra list comprehension:
In [ ]:
respuestas = [step_response(G, ts) for G in Gs]
y graficar, tambien, con una list comprehension:
In [ ]:
fig = figure(figsize=(12,3))
ax = fig.gca()
ps = [ax.plot(resp[0], resp[1]) for resp in respuestas]
ax.set_ylim(-0.01, 0.16)
ax.legend([ ps[0][0], ps[1][0], ps[2][0], ps[3][0], ps[4][0]],
["Estable", "Estable", "Crit. estable", "Inestable", "Inestable"]);
Si ahora tomamos uno de los sistemas inestables, podemos proponer un controlador que estabilice el sistema:
In [ ]:
ts = linspace(0, 10, 1000)
t, y = step_response(Gs[3], ts)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y)
por ejemplo, podemos proponer el controlador
$$ C(s) = K_P + s K_D $$en donde $K_P = 1$ y $K_D = 1$ como valores iniciales:
In [ ]:
ts = linspace(0, 10, 1000)
Con = tf([1, 1], [1])
t, y = step_response((Con*Gs[3]).feedback(), ts)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y)
Este controlador PD funciona para estabilizar el sistema, sin embargo aun no tiene las propiedades que deseo, por lo que tengo que sintonizar el sistema, para este proceso necesito cambiar $K_P$ y $K_D$ hasta obtener un comportamiento adecuado, propondré una notación alternativa para ayudarnos a escribir una ley de control mas compleja, primero analicemos el siguiente ejemplo:
$$ \frac{1}{s} + \frac{1}{s+1} = \frac{(s+1)}{s(s+1)} + \frac{s}{(s+1)s} = \frac{(s+1)+s}{s(s+1)} = \frac{2s+1}{s^2+s} $$
In [ ]:
tf([1], [1,0]) + tf([1], [1,1])
por lo que podemos concluir que si hacemos operaciones con las funciones de transferencia, el motor de algebra de las funciones de transferencia simplificará hasta el mínimo cualquier expresión, esto es resultado del hecho que las operaciones en Python, estan definidas como metodos de la clase que las define, por ejemplo, podemos preguntarle a un numero entero como se hace la suma en su clase:
In [ ]:
a = 1
In [ ]:
a.__add__??
o bien, a alguna de las funciones de transferencia:
In [ ]:
G1.__add__??
podriamos pensar hasta en sumar gatitos, pero esto quedará para otra práctica... por ahora conformemonos con definir el simbolo $s$, de tal manera que:
In [ ]:
s = tf([1, 0], [1])
s
In [ ]:
4*s + 1/s
Ahora definiremos nuestro controlador PID, el cual tiene la siguiente formula:
$$ C(s) = K_P + s K_D + \frac{1}{s} K_I $$
In [ ]:
kp, kd, ki = 10, 2, 1
Con = kp + kd*s + ki/s
Con
Pero como un solo ejemplo de controlador es muy aburrido, definiremos una función que simule el comportamiento de cualquier controlador PID que le demos:
In [ ]:
def sistema_controlado(kp, ki, kd):
Con = kp + kd*s + ki/s
ts = linspace(0, 10, 1000)
t, y = step_response((Con*Gs[3]).feedback(), ts)
fig = figure(figsize=(12,3))
ax = fig.gca()
ax.plot(t, y)
In [ ]:
sistema_controlado(1,1,1)
y recordando que podemos interactuar con estas variables, podemos hacer modificar nuestro controlador hasta que obtengamos un comportamiento adecuado:
In [ ]:
from ipywidgets import interact
In [ ]:
interact(sistema_controlado, kp=(0,20), ki=(0,20), kd=(0,20))
Sin embargo este sistema es lineal, y nuestros manipuladores robóticos no lo son; necesitaremos mas herramientas para simular el comportamiento de los robots, lo cual se verá en el documento llamado numerico.ipynb.