Modelizado Experimental¶
Vamos a encontrar el modelo matemático del siguiente sistema usando datos experimentales.
La planta es un sistema de iluminación con un LED de potencia (actuador) y un LDR (sensor). Este sistema es controlado a través de un Arduino Mega (actualmente Arduino Mega, en la foto Arduino Uno).
Toma de datos¶
Para la toma de datos usaremos el código en lazo abierto para Arduino. Y seguiremos los pasos mostrado a continuación.
- Descargamos el código.
- Subimos el código al Arduino.
- Abrimos el monitor serial.
- Copiamos los datos de un ciclo completo de apagado y prendido.
- Pegamos en un archivo llamado
datos.csvcon nombre de las columnast,u,y. - Finalmente, importamos los datos en el Colaboratory o Jupyter Notebook con pandas.
myHtml = """
<div>A continuación se muestran algunos datos. </div>
<link href="/assets/css/jupyter_nb.css" rel="stylesheet"/>
<script src="/assets/js/jupyter_nb.js"></script>
"""
display_html(myHtml, raw=True)
datos = pd.read_csv("datos.csv")
datos
| t | u | y | |
|---|---|---|---|
| 0 | 60968004 | 0 | 455 |
| 1 | 60972004 | 0 | 454 |
| 2 | 60976004 | 0 | 455 |
| 3 | 60980004 | 0 | 455 |
| 4 | 60984004 | 0 | 455 |
| ... | ... | ... | ... |
| 253 | 61980004 | 255 | 653 |
| 254 | 61984004 | 255 | 653 |
| 255 | 61988004 | 255 | 653 |
| 256 | 61992004 | 255 | 653 |
| 257 | 61996004 | 255 | 653 |
258 rows × 3 columns
Modelizado experimental¶
Despues de importar los datos, los extraemos en vectores y graficamos.
traw = np.array(datos["t"].tolist())
uraw = np.array(datos["u"].tolist())
yraw = np.array(datos["y"].tolist())
plt.plot(traw,uraw,traw,yraw);
plt.xlabel("tiempo (ms)");
plt.ylabel("bits");
plt.legend(["entrada","salida"]);
Para la construcción del modelo necesitamos tratar los datos, cambiaremos las unidades del tiempo que esta en microsegundos.
t = (traw - traw[0])/1e6
u = (uraw - uraw[0])
y = (yraw - yraw[0])
plt.plot(t,u,t,y);
plt.xlim((0,0.2));
plt.xlabel("tiempo (s)");
plt.ylabel("bits");
plt.legend(["entrada","salida"]);
Por la respuesta del sistema podemos suponer que es un sistema de primer orden con retardo. La función de transferencia para dicho sistema está dada por la expresión:
$$G(s) = \frac{\gamma}{\tau\,s+1}e^{-\theta\,s}$$Podemos suponer los valores de los parámetros y modificarlos hasta que las gráficas se parezcan.
s = tf("s")
gamma = y[-1]/u[-1]
tau = 0.0044376
theta = 0.004
num,den = pade(theta,n=5)
retardo = tf(num,den)
G = gamma/(tau*s+1)
print("------------------------")
print("Función de transferencia\n")
display(G)
GR = G*retardo
------------------------ Función de transferencia
La anterior expresión toma valores encontrados al tanteo.
Aproximación del modelo
T = np.linspace(t[0],t[-1],num=len(t))
ysim,tsim,_ = lsim(GR,U=u,T=T)
plt.plot(t,u,t,y,tsim,ysim);
plt.xlim((0,0.2));
plt.xlabel("tiempo (s)");
plt.ylabel("bits");
plt.legend(["entrada","salida","estimación"]);
Se ve que siguen habiendo diferencias entre el modelo y los datos. Aquí, podemos contruir un algoritmo de optimización que minimize el error cuadratico entre las dos respuestas (modelo y datos).
Estos son los resultados:
def model(x):
gamma, tau, theta = x
theta = max(abs(theta),0.001)
G = gamma/(tau*s+1)
num,den = pade(theta,n=5)
retardo = tf(num,den)
GR = G*retardo
ysim,tsim,_ = lsim(GR,U=u,T=T)
return np.sum((y-ysim)**2)
result = op.minimize(model,[0.77,0.005,0.004])
result
message: Desired error not necessarily achieved due to precision loss.
success: False
status: 2
fun: 171.0746620929305
x: [ 7.766e-01 5.023e-03 3.013e-03]
nit: 11
jac: [ 2.429e-01 1.100e+01 1.519e+01]
hess_inv: [[ 3.090e-08 2.338e-10 2.370e-10]
[ 2.338e-10 4.019e-11 4.061e-11]
[ 2.370e-10 4.061e-11 7.413e-11]]
nfev: 159
njev: 35
En el anterior resultado encontramos los valores optimizados para los parametros del modelo $\gamma$, $\tau$ y $\theta$. Con estos parámetros podemos verificar la respuesta nuevamente.
gamma, tau, theta = result.x
num,den = pade(theta,n=5)
retardo = tf(num,den)
G = gamma/(tau*s+1)
GR = G*retardo
print("Retardo:",theta)
Retardo: 0.003013355936680472
T = np.linspace(t[0],t[-1],num=len(t))
ysim,tsim,_ = lsim(GR,U=u,T=T)
plt.plot(t,u,t,y,tsim,ysim);
plt.xlim((0,0.2));
plt.xlabel("tiempo (s)");
plt.ylabel("bits");
plt.legend(["entrada","salida","estimación"]);
Diseño de un controlador (PID)¶
Para el diseño del controlador podemos usar un controlador PID o diseñar un controlador usando la representación en espacio de estados. Primero usaremos un PID, en este caso usamos un controlador PI, y los parámetros fueron ajustados manualmente.
Kp = 2.5
Ki = 3
Kd = 0.00
K = Kp + Ki/s + Kd*s
print("Controlador PI :")
display(K)
CPID = feedback(K*GR,1)
yy,tt = step(500*CPID)
plt.plot(tt,yy);
plt.xlabel("tiempo (s)");
plt.ylabel("bits");
Controlador PI :
print("Información de la respuesta a un escalón:\n")
stepinfo(500*CPID)
Información de la respuesta a un escalón:
{'RiseTime': 0.0034816640690269536,
'SettlingTime': 3.566964838718114,
'SettlingMin': 254.835424623772,
'SettlingMax': 500.74298048728696,
'Overshoot': 0.14859609745739136,
'Undershoot': 0,
'Peak': 500.74298048728696,
'PeakTime': 0.006963328138053907,
'SteadyStateValue': 500.0}
Se puede intentar optimizar el resultado utilizando una heuristica:
def controller(x):
Kp, Ki, Kd = x
Kp = max(0,Kp)
Ki = max(0,Ki)
Kd = max(0,Kd)
K = Kp + Ki/s + Kd*s
CPID = feedback(K*GR,1)
info = stepinfo(CPID)
return (info['Overshoot']+0.1)*(info['SettlingTime'])
result = op.minimize(controller,[2.5,3,0])
result
message: Desired error not necessarily achieved due to precision loss.
success: False
status: 2
fun: 0.3559594000412233
x: [ 2.490e+00 2.997e+00 -6.808e-03]
nit: 2
jac: [ 9.464e-02 2.333e+04 0.000e+00]
hess_inv: [[ 2.276e-01 3.415e-01 -4.008e-01]
[ 3.415e-01 5.122e-01 -6.012e-01]
[-4.008e-01 -6.012e-01 7.055e-01]]
nfev: 151
njev: 35
Con este resultado:
Kp, Ki, Kd = result.x
K = Kp + Ki/s
print("Controlador PI :")
display(K)
CPID = feedback(K*GR,1)
yy,tt = step(500*CPID)
plt.plot(tt,yy);
plt.xlabel("tiempo (s)");
plt.ylabel("bits");
Controlador PI :
