Modelizado Experimental¶
Encontrar modelos de sistemas a partir de datos experimentales y probar controladores en dichos sistemas.
 |
 |
 |
 |
 |
 |
 |
Modelizado experimental de un sistema de iluminación¶
Para entender el modelizado experimental usaremos una maqueta (planta) para el control de iluminación.
La planta es un sistema de iluminación con un LED Light Emitter Diode de potencia (actuador) y un LDR Light Dependent Resistor (sensor). Este sistema es controlado a través de un Arduino X.
Encontrar el modelo de forma analitica para este sistema puede ser dificil dado que no todos tenemos conocimientos en óptica. Por lo que para encontrar el modelo matemático de este, usaremos datos experimentales.
1. Subir el Código Arduino¶
Abramos el archivo iluminacion.ino) en el Arduino IDE (IDE Integrated Development Environment). Este código nos permitira actuar el LED y obtener los datos del sensor. Para interactuar con el arduino usaremos el puerto serial y python.
Para ubicar el puerto de comunicación serial vamos al menú Herramientas > Puerto. Seleccionamos el puerto donde se encuentre el Arduino X.
Verificamos en la maqueta donde se tiene conectada la entrada y la salida. Y corregimos los pines en el código.
#define OUT 2
#define IN A0
Luego, subimos el código al Arduino X.
Este código básicamente lee la información del sensor y la envía por puerto serial al computador, esta información esta codificada en código ASCII. Adicionalmente si uno envía un valor (0-255) por el puerto serial al Arduino, este cambiará la intensidad del LED proporcionalmente a este valor.
Una vez subido el código al Arduino se puede verificar vía el Monitor Serie la recepción de caracteres
ASCII, enviando cualquier carácter, el carácter recibido debería cambiar.
2. Análisis de la toma de datos¶
Despues de ejecutar el código de python que genera cambios en el actuador, tenemos la siguiente respuesta.
plt.plot(T,X)
plt.plot(T,Y)
plt.show()
Vamos a filtrar la respuesta con un promedio (average).
t = numpy.array(T)
x = numpy.array(X)
y = numpy.array(Y)
x = x.astype(float)
y = y.astype(float)
y = y - numpy.average(y[0:50])
plt.plot(t,x)
plt.plot(t,y)
plt.show()
Hagamos zoom en donde se genera el escalón:
plt.plot(t[90:110],x[90:110])
plt.plot(t[90:110],y[90:110])
plt.show()
Respodamos las siguientes preguntas:
- ¿El sistema es de primer o segundo orden?
- ¿Cuál es su ganancia $\gamma$?
- ¿Cuál es su frecuencia angular natural $\omega_n$?
- ¿Cuál es su coefficiente de amortiguamiento $\xi$?
- ¿Cuál es su constante de tiempo?
3. Indentificación paramétrica¶
Para encontrar los paramétros del sistema.
$$\frac{\gamma}{\tau\, s + 1}\, e^{-\theta\,s}$$debemos referirnos a la curva experimental.
- La ganancia puede ser calculada como: $$\gamma = \frac{y_\infty - y_0}{u_\infty - u_0}$$
- El retardo $\theta$ es el tiempo que se demora la salida en responder ante un cambio de la entrada en el sistema.
- El tiempo de establecimieto $\tau$ es el tiempo que se demora la respuesta en alcazar un $63.2\%$ del su valor final.
Con los parámetros obtenidos podemos simular el sistema y verificar que tan cerca esta del comportamiento real.
identificacion(0.76,0.003,0.005)
plt.xlim(0.95,1.05);
4. Optimización paramétrica¶
Por otro lado los parametros que encontramos manualmente, se puede buscar utilizando optimizadores que ya estan disponibles en python, en la libreria Scipy.
Optimizando los paramátros del nuestro sistema experimental tenemos:
def simulador(x,gamma,theta,tau):
fo = control.tf([gamma],[tau,1])
dt_num,dt_den = control.pade(theta,1)
delay = control.tf(dt_num,dt_den)
fodt = fo*delay
tin = numpy.linspace(0.01, 3, len(t))
tsim, ysim, xsim = control.forced_response(fodt, tin, x, X0= 0.0)
for i in range(0,len(ysim)):
ysim[i] = max(ysim[i],0)
return ysim
popt, pcov = scipy.optimize.curve_fit(simulador, x, y, bounds=([0.001,0.001,0.001], [3., 1., 0.5]))
gamma = popt[0]
theta = popt[1]
tau = popt[2]
print(" - gamma = %2.4f\n - tau = %2.4f\n - theta = %2.4f" %(gamma,tau,theta))
ysim = simulador(x,gamma,theta,tau)
Dando como resultado:
resultado_optimizacion(0.95,1.05)
5. Análisis de estabilidad del resultado¶
Habiendo identificado el sistema, podemos analizar su estabilidad. Para esto podemos hacer uso del root_locus.
fo = control.tf([gamma],[tau,1])
dt_num,dt_den = control.pade(theta,1)
delay = control.tf(dt_num,dt_den)
fodt = fo*delay
control.root_locus(fodt);
plt.axvspan(0, 2000, facecolor='r', alpha=0.5)
plt.xlim(-2000,2000);
Como nos muestra la gráfica, existen valores de $K$ para los cuales el sistema es inestable.
6. Propuesta de controlador¶
Como controlador podemos usar un control PID o podemos realizar diseñar un controlador por posicionamiento de polos.
Empecemos con un control PID, para este vamos a encontrar la ganancia ultima del sistema $K_u$ y el periodo ultimo $T_u$
import math
pi = math.pi
Ku,_,Fu,_ = control.margin(fodt)
Tu = 2*pi/Fu
print("La ganancia ultima para el sistema experimental es Ku = %2.4f y el periodo último es Tu = %2.4f"%(Ku,Tu))
Con estos valores podemos usar la tabla de sintonia de Ziegler-Nichols. Para un controlador PI tenemos:
s = control.tf([1,0],1)
kp = 0.45*Ku
ti = Tu/1.2
td = 0
K = kp * (1 + 1/(ti*s)+td*s)
display(K)
Podemos verificar la estabilidad del sistema con este controlador, para esto encontremos los polos del sistema en lazo cerrado.
G = K * fodt
closed_loop = control.feedback(G,1)
polos = control.pole(closed_loop)
print("Los polos del sistema son: \n\n%s" % ('\n'.join(map(str, polos))))
Como se puede observar todos tiene la parte real negativa, por lo tanto es sistema es estable.
7. Implementación del controlador¶
Para implementar el PID en un sistema de tiempo discreto como lo es un Arduino u otro microcontrolador. Debemos discretizar el controlador. Esto se hace con la transformada $z$, remplazando:
$$s \to \frac{2}{T}\frac{z-1}{z+1}$$donde $T$ es el tiempo de muestreo.
En python este trabajo lo logramos con sample_system:
Kdiscreto = control.sample_system(K,0.01)
display(Kdiscreto)
De esta función de transferencia discreta debemos sacar el código que usaremos en el Arduino. Para esto debemos construir una ecuación que tenga todos los terminos con $z^n$, con $n$ negativo o zero, de ahí, despejamos el valor del actuador que esta acompañado por $z^0$
z,U,UF,E = sympy.symbols('z U U_k E')
numer = Kdiscreto.num[0][0]
denom = Kdiscreto.den[0][0]
rhs = 0
for i,d in enumerate(denom):
if i == 0:
rhs += UF*d*z**(len(denom)-i-1)
else:
rhs += U*d*z**(len(denom)-i-1)
lhs = 0
for i,n in enumerate(numer):
lhs += E*n*z**(len(numer)-i-1)
tfd = sympy.Eq(sympy.expand(rhs/z**(len(numer)-1)),sympy.expand(lhs/z**(len(numer)-1)))
codigo = sympy.Eq(UF,sympy.expand(sympy.solve(tfd,UF)[0]))
display(codigo)
La anterior expresión es la expresión que se programa para el controlador.
// CONTROLADOR
E0 = setpoint - float(sensor) // Calculo del Error
U0 = 0.72*E0 + 0.13*E1 + 1*U1 // Calculo de la expresion para el actuador
E1 = E0 // Guardo el error actual en la variable de error anterior
U1 = U0 // Hago lo mismo para el actuador
8. Verificacion del controlador¶
Respuesta del sistema al controlador propuesto
def show_results():
t = numpy.array(TC)
x = numpy.array(XC)
y = numpy.array(YC)
s = numpy.array(SC)
x = x.astype(float)
y = y.astype(float)
s = s.astype(float)
ax = plt.subplot(111)
ax.plot(t[1:],s[1:], label="respuesta")
ax.plot(t[1:],y[1:], label="referencia")
#ax.plot(t[1:],x[1:], label="actuador")
ax.legend()
plt.title('Respuesta del sistema en lazo cerrado')
plt.xlabel('tiempo (t)')
plt.ylabel('respuesta')
show_results()
Otro ejemplo de respuesta
show_results()
