Control de un reactor continuo de tanque agitado
Abstract
En este trabajo se desarrolla una aplicación informática con una interfaz gráfica de usuario que facilitará el aprendizaje de los fundamentos del control automático aplicado a procesos químicos, específicamente, en reactores de tanque agitado con dinámica no lineal. El desarrollo se ha hecho mediante la programación en Matlab como entorno de simulación. Para ello, en primer lugar, se han obtenido tres modelos dinámicos de un reactor de tanque agitado continuo con diferentes métodos de refrigeración: mediante un serpentín y mediante una camisa con y sin recirculación de refrigerante. Posteriormente, se ha programado un simulador para integrar las ecuaciones de los modelos a partir de unas condiciones iniciales, mediante los métodos de Euler y Runge-Kutta de orden 4. Asimismo, se ha usado un temporizador para permitir la simulación en tiempo real o en tiempo de simulación, con una velocidad del proceso más rápida. Por otro lado, se ha diseñado un algoritmo para la identificación de la dinámica del proceso. El ajuste de los parámetros se realiza aplicando el método de mínimos cuadrados, sobre modelos lineales de primer orden con tiempo muerto, partiendo de los datos obtenidos por simulación. Estos modelos simplifican la dinámica no lineal del proceso real. Para cada modelo, se han programado las ecuaciones de distintos reguladores, con el objetivo de controlar la temperatura del reactor. Se han utilizado algoritmos PID con y sin acción anticipativa, realimentación del estado y control predictivo DMC monovariable. Además, para el reactor encamisado, se han desarrollado las estructuras de control en cascada, control adaptativo de ajuste por tabla y control autoajustable basado en mínimos cuadrados recursivos. Todos los parámetros de ajuste han sido diseñados para proporcionar al usuario una solución por defecto. Finalmente, se han integrado todas estas funcionalidades en una interfaz gráfica de usuario multiventana, que permite el ajuste de las diferentes opciones de simulación y la selección del algoritmo de control. Además, se ha incluido una ventana para realizar ensayos de entrada en escalón, permitiendo exportar a una hoja de cálculo los datos obtenidos para realizar la identificación. Con fines didácticos, se ha creado un manual en el que se detalla el procedimiento a seguir para la instalación y manejo del software diseñado.
Full text
Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Control de un reactor continuo de tanque agitado Autor: Jesús Fernández Romero Tutor: Francisco Javier Gutiérrez Ortiz Dpto. de Ingeniería Química y Ambiental Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021
i Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Control de un reactor continuo de tanque agitado Autor: Jesús Fernández Romero Tutor: Francisco Javier Gutiérrez Ortiz Catedrático de Universidad Dpto. de Ingeniería Química y Ambiental Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022
iii Trabajo Fin de Máster: Control de un reactor continuo de tanque agitado Autor: Jesús Fernández Romero Tutor: Francisco Javier Gutiérrez Ortiz El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2022 El Secretario del Tribunal
v Agradecimientos Quiero agradecer este Trabajo Fin de Máster a: Mis profesores, por su trabajo y dedicación. Mi tutor, por su constante esfuerzo dirigiéndome este TFM, respondiéndome incluso los fines de semana. Mis amigos, por ayudarme en los malos momentos que he tenido. Mi familia, por ser un pilar fundamental en mi vida, apoyándome siempre en todas mis decisiones. En especial, quiero dedicarle “In memoriam” este Trabajo Fin de Master a mi hermano Alejandro, quien siempre ha estado y estará a mi lado en todo momento. Jesús Fernández Romero Sevilla, 2022
vii Resumen En este trabajo se desarrolla una aplicación informática con una interfaz gráfica de usuario que facilitará el aprendizaje de los fundamentos del control automático aplicado a procesos químicos, específicamente, en reactores de tanque agitado con dinámica no lineal. El desarrollo se ha hecho mediante la programación en Matlab como entorno de simulación. Para ello, en primer lugar, se han obtenido tres modelos dinámicos de un reactor de tanque agitado continuo con diferentes métodos de refrigeración: mediante un serpentín y mediante una camisa con y sin recirculación de refrigerante. Posteriormente, se ha programado un simulador para integrar las ecuaciones de los modelos a partir de unas condiciones iniciales, mediante los métodos de Euler y Runge-Kutta de orden 4. Asimismo, se ha usado un temporizador para permitir la simulación en tiempo real o en tiempo de simulación, con una velocidad del proceso más rápida. Por otro lado, se ha diseñado un algoritmo para la identificación de la dinámica del proceso. El ajuste de los parámetros se realiza aplicando el método de mínimos cuadrados, sobre modelos lineales de primer orden con tiempo muerto, partiendo de los datos obtenidos por simulación. Estos modelos simplifican la dinámica no lineal del proceso real. Para cada modelo, se han programado las ecuaciones de distintos reguladores, con el objetivo de controlar la temperatura del reactor. Se han utilizado algoritmos PID con y sin acción anticipativa, realimentación del estado y control predictivo DMC monovariable. Además, para el reactor encamisado, se han desarrollado las estructuras de control en cascada, control adaptativo de ajuste por tabla y control autoajustable basado en mínimos cuadrados recursivos. Todos los parámetros de ajuste han sido diseñados para proporcionar al usuario una solución por defecto. Finalmente, se han integrado todas estas funcionalidades en una interfaz gráfica de usuario multiventana, que permite el ajuste de las diferentes opciones de simulación y la selección del algoritmo de control. Además, se ha incluido una ventana para realizar ensayos de entrada en escalón, permitiendo exportar a una hoja de cálculo los datos obtenidos para realizar la identificación. Con fines didácticos, se ha creado un manual en el que se detalla el procedimiento a seguir para la instalación y manejo del software diseñado.
xv Índice de Figuras Figura 1. Esquema de un reactor de tanque agitado continuo y refrigerado mediante un serpentín ................ 3 Figura 2. Esquema de un reactor de tanque agitado continuo con camisa de refrigeración .............................. 7 Figura 3. Esquema de un reactor de tanque agitado continuo encamisado con recirculación de refrigerante . 9 Figura 4. Diagrama de flujo del algoritmo de simulación dinámica en Matlab ................................................. 19 Figura 5. Diagrama de bloques de un modelo con estructura ARX ................................................................... 22 Figura 6. Diagrama de bloques de un modelo con estructura ARMAX ............................................................. 23 Figura 7. Diagrama de bloques de un modelo con estructura Box and Jenkins ................................................ 24 Figura 8. Suavizado de datos mediante filtro de media móvil ........................................................................... 26 Figura 9. Comparación entre proceso modelado y proceso identificado empíricamente para distintos modelos de reactor CSTR ........................................................................................................................... 27 Figura 10. Diagrama de bloques del sistema lineal completo ........................................................................... 29 Figura 11. Realización en la GUI de una identificación paramétrica por cada punto de operación del ensayo 31 Figura 12. Herramienta de análisis y modificación de funciones de transferencia ........................................... 33 Figura 13. Comparación de desempeño entre distintos tipos de acción derivativa .......................................... 39 Figura 14. Diagrama de bloques de la implementación en la GUI implementado en la GUI ............................ 40 Figura 15. Comparación de desempeño entre distintos controladores PI obtenidos mediante sintonización IMC ............................................................................................................................................................. 42 Figura 16. Diagrama de bloques del control en lazo cerrado con acción anticipativa ...................................... 43 Figura 17. Comparación de las respuestas temporales ante una perturbación para los diferentes modelos controlados por un PI con y sin acción anticipativa ................................................................................... 44 Figura 18. Implementación en la GUI del control en lazo cerrado con acción anticipativa ............................... 45 Figura 19. Diagrama de bloques del algoritmo de control en cascada .............................................................. 46 Figura 20. Rechazo de una perturbación (temperatura del refrigerante a la entrada de la camisa del reactor) utilizando diferentes variaciones del algoritmo PI..................................................................................... 47 Figura 21. Rechazo de una perturbación externa (temperatura a la entrada del reactor) para distintas estrategias de control ................................................................................................................................ 47 Figura 22. Desempeño de un controlador “state feedback” proporcional para distintos valores de 𝜇 ............ 51 Figura 23. Desempeño de un controlador “state feedback” con acción integral para distintos valores de 𝛾 .. 51 Figura 24. Desempeño en el rechazo de perturbaciones de un controlador “state feedback” con acción integral para distintos valores de 𝛾 ........................................................................................................... 52 Figura 25. Herramienta para la sintonización de controladores de realimentación de estado ........................ 53 Figura 26. Diagrama de funcionamiento de un controlador predictivo basado en modelo ............................. 55 Figura 27. Desempeño de un controlador DMC para distintos valores de 𝜆 .................................................... 60
xvi Figura 28. Comparación del rechazo de una perturbación para distintas variaciones del algoritmo DMC ....... 60 Figura 29. Comparación de la respuesta temporal obtenida mediante control PID y DMC .............................. 61 Figura 30. Diagrama de bloques de un control adaptativo mediante ajuste por tabla ..................................... 64 Figura 31. Diagrama de conmutación de la estrategia de control adaptativo de ajuste por tabla .................... 65 Figura 32. Comparación entre las respuestas temporales obtenidas para el caso del reactor CSTR encamisado, aplicando control PI y DMC con y sin estrategia de ajuste por tabla. .................................. 67 Figura 33. Diagrama de bloques de un controlador autoajustable ................................................................... 68 Figura 34. Resultados de una identificación paramétrica mediante mínimos cuadrados recursivos ............... 70 Figura 35. Respuesta temporal y evolución de los parámetros 𝑘𝑐 y 𝑇𝑖 para un control PI autoajustable, aplicado al caso del reactor CSTR encamisado .......................................................................................... 73 Figura A1. Instalación del software. Inicio.......................................................................................................... 79 Figura A2. Instalación del software. Selección de ruta de instalación ............................................................... 80 Figura A3. Instalación del software Matlab Runtime ......................................................................................... 80 Figura A4. Instalación del software. Resumen de la instalación ........................................................................ 81 Figura A5. Ventana de inicio de la interfaz gráfica de usuario ........................................................................... 81 Figura A6. Ventana principal de la interfaz gráfica de usuario ........................................................................... 82 Figura A7. Evolución del estado interno del proceso obtenido mediante el botón “Variables de proceso” .... 83 Figura A8. Evolución del error y de la integral del error obtenidos mediante el botón “Índices de error” ....... 83 Figura A9. Evolución de las señales eléctricas obtenidas mediante el botón “Señales eléctricas” ................... 84 Figura A10. Ventana gráfica de configuración de la simulación ........................................................................ 85 Figura A11. Opciones de gráficas ....................................................................................................................... 87 Figura A12. Guardado de datos en el modo de simulación dinámica ............................................................... 88 Figura A13. Ventana gráfica de identificación paramétrica ............................................................................... 89 Figura A14. Importación de datos ...................................................................................................................... 90 Figura A15. Mensaje de confirmación de exportación correcta de datos y directorio de ubicación ................ 91 Figura A16. Resultados de la identificación paramétrica mediante la opción “Identificación por puntos” ...... 92 Figura A17. Modelos lineales obtenidos con una identificación paramétrica mediante la opción “Identificación por puntos” ........................................................................................................................ 93 Figura A18. Herramienta para el análisis de funciones de transferencia .......................................................... 94 Figura A19. Ventana gráfica de opciones de control ......................................................................................... 96 Figura A20. Previsualización de la respuesta temporal del proceso controlado para una referencia unitaria . 97 Figura A21. Ventana gráfica de información sobre el controlador seleccionado .............................................. 97 Figura A22. Sintonización de un controlador PID............................................................................................... 98 Figura A23. Configuración del tipo de acción proporcional ............................................................................... 99 Figura A24. Configuración del tipo de acción derivativa .................................................................................... 99 Figura A25. Configuración del tipo de acción derivativa con filtro .................................................................. 100 Figura A26. Diagrama de bode del proceso controlado en lazo abierto (I) ..................................................... 100 Figura A27. Diagrama de bode del proceso controlado en lazo abierto (II) .................................................... 101 Figura A28. Sintonización de un controlador PI con acción anticipativa ......................................................... 102
xvii Figura A29. Herramienta para sintonizar la acción anticipativa ...................................................................... 103 Figura A30. Lazo de control configurado con acción anticipativa ................................................................... 104 Figura A31. Sintonización de un controlador por realimentación del estado con acción integral .................. 104 Figura A32. Herramienta para sintonización de un controlador por realimentación del estado.................... 105 Figura A33. Sintonización de un controlador predictivo DMC ........................................................................ 106 Figura A34. Ventana gráfica de opciones de control para el reactor encamisado con estructura de control en cascada..................................................................................................................................................... 107 Figura A35. Ventana gráfica de control adaptativo y autoajustable ............................................................... 109 Figura A36. Parámetros de ajuste para el control PI con estrategia “Gain Scheduling” ................................. 110 Figura A37. Parámetros de ajuste para el control autoajustable basado en RLS ............................................ 111 Figura A38. Parámetros de ajuste para el control predictivo DMC con estrategia “Gain Scheduling” ........... 111 Figura B1. Resultados de la identificación paramétrica LS .............................................................................. 115 Figura B2. Resultados de la identificación paramétrica RLS y comparación con el método LS ....................... 119
xix Notación Símbolos 𝑎 Una barra localizada sobre una letra indica un vector 𝑏 Un punto sobre una letra indica una variable derivada con respecto al tiempo 𝐶𝑒𝑖𝑙(𝑥) Redondeo al entero más próximo por exceso 𝑐𝑡𝑒 Valor constante 𝑓(𝑥) |𝑥=𝑥0 Evaluación en el punto de equilibrio de una función de 𝑥 𝐼 Matriz identidad 𝐽∗(𝑥) Valor óptimo de la función de coste 𝐽(𝑥) 𝐽∗(𝑥+) Valor óptimo de la función de coste 𝐽(𝑥), analizado en 𝑥(𝑘+1). ln(·) Logaritmo neperiano 𝑀 Una doble barra localizada sobre una letra mayúscula indica una matriz 𝑀𝑗,∗ Fila 𝑗 de la matriz 𝑀 𝑀(𝑖,𝑗:𝑘) Elementos de M localizados en la fila 𝑖, entre las columnas desde 𝑗 hasta 𝑘 min𝐽(𝑥) Valor mínimo de la función 𝐽(𝑥) 𝑃0 Valor inicial de la variable 𝑃 𝑠 Variable compleja usada en la transformada de Laplace 𝑥∗ Valor que optimiza a una función objetivo 𝐽(𝑥) 𝑧 Variable compleja usada en la transformada 𝒵 𝒵{·} Transformada 𝒵 de la función {·} Δ Incremento (·)𝑇 Transpuesta de (·) |·| Valor absoluto de (·) 𝜉(·) Ruido blanco de media cero y varianza 𝜎2 𝜙(𝑠) Polinomio característico evaluado en 𝑠
xx Variables y propiedades fisicoquímicas 𝐴 Área del reactor (𝑚2) 𝑐𝐴 Concentración molar de reactivo (𝑚𝑜𝑙/𝑚3) 𝐶𝑝 Capacidad calorífica a presión constante del reactivo y producto (𝑘𝐽/(𝑘𝑔 ℃) ) 𝐸 Energía de activación (𝐽/𝑚𝑜𝑙) 𝐹 Caudal de reactivo y producto (𝑚3/𝑠) 𝑘 Constante cinética (𝑠−1) 𝑘0 Factor de frecuencia (𝑠−1) 𝐾𝑣𝑠 Coeficiente de caudal máximo de la válvula (𝑚3ℎ−1𝑏𝑎𝑟−1/2) 𝑄 Potencia calorífica liberada hacia el refrigerante (𝑊) 𝑟 Velocidad de reacción (𝑚𝑜𝑙/(𝑚3 𝑠)) 𝑅 Constante universal de los gases ideales (8,314 𝐽/(𝑚𝑜𝑙 𝐾)) 𝑇 Temperatura del reactor (℃) 𝑈 Coeficiente global de transferencia de calor (𝑊/(𝑚2 ℃)) 𝑉 Volumen del reactor (𝑚3) 𝑉𝑐 Volumen de la camisa (𝑚3) Δ𝐻𝑅 Calor de reacción (𝐽/𝑚𝑜𝑙) Δ𝑃𝑣(𝑥) Caída de presión en la válvula (𝑏𝑎𝑟) 𝛾 Peso específico del refrigerante 𝜌 Densidad del reactivo y producto (consideradas iguales) (𝑘𝑔/𝑚3)
xxi Variables de modelos lineales e identificación 𝐴 Matriz de estados del modelo en espacio de estados 𝐵 Matriz de entrada del modelo en espacio de estados 𝐶 Matriz de salida del modelo en espacio de estados 𝐷 Matriz de transmisión directa del modelo en espacio de estados 𝐸 Vector de error de modelado 𝑓𝑝𝑐𝑡 Sensibilidad del bloque conversor de la salida del transmisor (en 𝑚𝐴 o 𝑉) a un valor porcentual 𝑓𝑐 Conversión de la señal de control calculada en señal eléctrica (de % a 𝑉 o 𝑚𝐴) 𝑓𝑠𝑡 Conversión de la señal eléctrica en magnitud física 𝑓 Respuesta libre 𝑔 Vector de ganancias del modelo de respuesta en escalón 𝑔𝑑 Vector de ganancias del modelo de respuesta en escalón para una entrada de perturbación 𝐺 Matriz dinámica 𝐺 Función de transferencia 𝐾 Ganancia estática 𝑀𝐶 Matriz de controlabilidad 𝑁1 Longitud del vector de estados 𝑛𝑎 Orden de la parte autorregresiva de una estructura de modelo lineal en tiempo discreto 𝑛𝑏 Orden de la parte exógena de una estructura de modelo lineal en tiempo discreto 𝑛𝑐 Orden de la parte de media móvil de una estructura de modelo lineal en tiempo discreto 𝑛𝑑𝑒𝑙 Número de pasos que tarda un cambio en la entrada en afectar a la señal de salida 𝑥 Vector de estados 𝑋 Base de datos de entrada al proceso 𝑦(𝑡+𝑘 | 𝑡) Valor predicho de la variable "𝑦" en el instante "𝑡+𝑘" calculado en el instante actual, "𝑡", con 𝑘∈ℕ+. Cada número natural se corresponde con un intervalo de tiempo igual al periodo de muestro 𝑦 Salida del proceso 𝑌 Base de datos de salida al proceso Γ Matriz de varianza-covarianza 𝜖 Error de modelado 𝜁 Constante de amortiguamiento 𝜃 Vector paramétrico 𝜆𝑜 Factor de olvido 𝜏 Constante de tiempo (𝑚𝑖𝑛) 𝜏𝑚 Tiempo muerto 𝜐(𝑡) Término de media móvil asociado al ruido en una estructura de tipo ARMAX 𝜑 Vector de regresión
xxii Variables de simulación y control 𝑑 Vector de perturbación 𝐷𝑦 Derivada de la salida del proceso con respeto al tiempo 𝑒𝑟 Error de seguimiento 𝑢𝐷 Término de la acción derivativa de una estructura PID estándar 𝐻 Matriz de históricos 𝐾𝑐 Matriz de ganancia de realimentación 𝐾𝐼 Elementos de la matriz de ganancia de realimentación correspondientes a la integral del error 𝐼𝐸 Integral del error 𝐾𝑐 Constante proporcional de una estructura PID estándar 𝐾𝑐𝑑 Constante proporcional de una estructura PID estándar en tiempo discreto 𝑚 Horizonte de control 𝑁2 Longitud del vector de referencia 𝑁3 Longitud del vector de perturbación 𝑛1 Número de filas de la matriz de históricos 𝑛2 Número de columnas de la matriz de históricos 𝑝 Horizonte de predicción 𝑃𝑟𝑒𝐶 Precompensador 𝑄 Matriz de ponderación de desviación respecto a la posición de equilibrio 𝑟 Vector de referencias 𝑅 Matriz de ponderación de la energía de control 𝑡𝑑 Vector de instantes de aplicación de las perturbaciones en la simulación 𝑡𝑟 Vector de instantes de aplicación de las referencias en la simulación 𝑇𝑑 Tiempo derivativo (𝑚𝑖𝑛) 𝑇𝑑𝑑 Tiempo derivativo de una estructura PID estándar en tiempo discreto (𝑚𝑖𝑛) 𝑇𝐷𝑀𝐶 Tiempo de muestreo del algoritmo de control DMC (𝑚𝑖𝑛) 𝑇𝑖 Tiempo integral (𝑚𝑖𝑛) 𝑇𝑖𝑑 Tiempo integral de una estructura PID estándar en tiempo discreto (𝑚𝑖𝑛) 𝑡𝑚 Tiempo de muestreo (𝑚𝑖𝑛) 𝑇𝑠 Tiempo de simulación (𝑚𝑖𝑛) 𝑢 Señal de control 𝑢𝑜𝑠 Valor nominal de salida de un controlador proporcional 𝑤 Vector de referencias futuras Δ𝑡𝑘 Paso de integración (min) 𝜆 Parámetro de ponderación de los esfuerzos de control
xxiii Subíndices 𝑎 El subíndice 𝑎 sobre una función de transferencia hace referencia al actuador (p.e. 𝐺𝑎𝑣) 𝑐 El subíndice 𝑐 sobre un parámetro o variable fisicoquímica hace referencia al caudal de refrigeración (p.e. 𝐹𝑐) 𝑑 El subíndice 𝑑 sobre una variable de modelo lineal hace referencia a la perturbación (p.e. 𝐺𝑑) 𝑑𝑒𝑠 El subíndice 𝑑𝑒𝑠 sobre un parámetro o variable hace referencia al valor deseado (p.e. 𝐾𝑐,𝑑𝑒𝑠) 𝑒 El subíndice 𝑒 sobre un parámetro o variable fisicoquímica hace referencia la entrada del sistema (p.e. 𝑇𝑐𝑒) 𝑗 El subíndice 𝑗 sobre un parámetro o variable fisicoquímica hace referencia al caudal de mezcla (p.e. 𝑇𝑗) 𝑘 El subíndice 𝑘 sobre una variable de simulación hace referencia al valor en el instante actual (p.e. 𝑟𝑘) 𝑚á𝑥 El subíndice 𝑚á𝑥 sobre una variable hace referencia al valor máximo alcanzable o valor máximo establecido (p.e. 𝑝𝑚á𝑥) 𝑜𝑠 El subíndice 𝑜𝑠 sobre una variable de control hace referencia al valor de “offset” 𝑜𝑢𝑡 El subíndice 𝑜𝑢𝑡 sobre una variable de simulación hace referencia al argumento de salida del algoritmo de simulación (p.e. 𝐻𝑜𝑢𝑡) 𝑝 El subíndice 𝑝 sobre una variable de modelo lineal hace referencia al proceso (p.e. 𝐺𝑝) 𝑟 El subíndice 𝑟 sobre un parámetro o variable fisicoquímica hace referencia al caudal recirculado (p.e. 𝐹𝑟𝑐) 𝑟𝑒𝑞 El subíndice 𝑟𝑒𝑞 sobre un parámetro o variable fisicoquímica hace referencia al caudal requerido (p.e. 𝐹𝑐,𝑟𝑒𝑞) 𝑠 El subíndice 𝑠 sobre un parámetro o variable fisicoquímica hace referencia al caudal de purga (p.e. 𝐹𝑐𝑠). El subíndice 𝑠 sobre una constante de tiempo hace referencia una función de transferencia de segundo orden. 𝑇 El subíndice 𝑇 sobre una función de transferencia hace referencia a la combinación de las dinámicas dadas por 𝑓𝑐, 𝐺𝑎𝑣, 𝐺𝑣 y 𝐺𝑝 (p.e. 𝐺𝑇) 𝑇𝑑 El subíndice 𝑇𝑑 sobre una función de transferencia hace referencia al transmisor (p.e. 𝐺𝑇𝑑) 𝑣 El subíndice 𝑣 sobre una variable de modelo lineal hace referencia a la válvula de control (p.e. 𝐺𝑣)
Modelado del reactor de tanque agitado continuo 4 Para obtener un modelo matemático sencillo es necesario realizar algunas simplificaciones, que se detallan a continuación: La pureza del reactivo a la entrada del reactor es del 100%. La extracción de calor de la reacción exotérmica se lleva a cabo mediante un serpentín por el que circula un caudal de líquido de refrigeración. La mezcla dentro del reactor está perfectamente agitada. La densidad del reactivo a la entrada es igual a la densidad del producto a la salida. Existe un control de nivel en el tanque asociado al caudal de salida, de forma que el volumen de la mezcla en el interior se mantiene constante en el tiempo. Se desprecian las pérdidas de calor hacia el exterior. Los calores específicos a presión constante (𝐶𝑝) no cambian con la temperatura dentro del intervalo de operación, y coinciden en las corrientes del proceso a la entrada y salida del reactor. De esta forma, es posible realizar un balance de masa en el reactor: 𝑉𝑑𝑐𝐴 𝑑𝑡=𝐹(𝑐𝐴𝑒−𝑐𝐴)−𝑉𝑘𝑐𝐴 (3) El calor cedido por el reactor hacia el refrigerante por unidad de tiempo (potencia calorífica) puede modelarse como: 𝑄=𝜌𝑐𝐶𝑝𝑐𝐹𝑐((𝑇𝑐−𝑇𝑐𝑒)=𝑈𝐴 𝑇𝑐−𝑇𝑐𝑒 𝐿𝑛(𝑇−𝑇𝑐𝑒 𝑇−𝑇𝑐) (4) de donde se obtiene: 𝐿𝑛(𝑇−𝑇𝑐𝑒 𝑇−𝑇𝑐)= 𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐 (5) Por lo que 𝑄 puede expresarse como: 𝑄=𝜌𝑐𝐶𝑝𝑐𝐹𝑐(1−𝑒−𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐)(𝑇−𝑇𝑐𝑒) (6) Realizando un balance de energía se obtiene la siguiente ecuación diferencial: 𝑉𝜌𝐶𝑝𝑑𝑇 𝑑𝑡=𝜌𝐶𝑝𝐹(𝑇𝑒−𝑇)+(−Δ𝐻𝑅)𝑉𝑘𝑐𝐴−𝜌𝑐𝐶𝑝𝑐𝐹𝑐(1−𝑒−𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐)(𝑇−𝑇𝑐𝑒) (7) Las ecuaciones (3) y (7) componen un sistema de dos ecuaciones diferenciales que describen la dinámica del reactor CSTR con refrigeración por serpentín. Las variables de salida son la temperatura del reactor y la concentración del reactivo A en el caudal saliente, y las entradas son el caudal de refrigerante, su temperatura de entrada, la temperatura de la corriente de alimentación al reactor, el caudal del proceso y la concentración a la entrada. En la Tabla 1 se muestran los parámetros del modelo implementado en la interfaz gráfica, así como las condiciones iniciales. La temperatura del refrigerante y de la corriente de alimentación, así como el caudal del proceso y la concentración de reactivo a la entrada, son posibles perturbaciones.
5 Control de un reactor continuo de tanque agitado Considerando que el caudal de refrigeración puede manipularse, dentro de un determinado intervalo, mediante una válvula de control automática situada a la entrada del serpentín, se tiene: 𝐹𝑐=𝑓(𝑥)𝐾𝑣𝑠√Δ𝑃𝑣(𝑥) 𝛾 (8) donde: 𝐾𝑣𝑠 es el coeficiente de caudal máximo de la válvula y se obtiene cuando 𝑥=1 y 𝑓(𝑥)=1. 𝑓(𝑥) es la característica inherente de la válvula. 𝐾𝑣=𝐾𝑣𝑠𝑓(𝑥) es el coeficiente de válvula; varía con la abertura de la válvula en el intervalo 𝑥∈[0,1]. Δ𝑃𝑣(𝑥) es la caída de presión en la válvula. 𝛾 es el peso específico del refrigerante (se tomará igual a 1, por tratarse de agua). Realizando las siguientes simplificaciones: Δ𝑃𝑣(𝑥)=𝑐𝑡𝑒 (9) 𝑓(𝑥)=𝑥 (10) La ecuación de la válvula puede aproximarse como sigue: 𝐹𝑐=𝑥·𝐹𝑐,𝑚á𝑥 (11) De forma que puede establecerse una relación entre la abertura de la válvula y el caudal solicitado correspondiente a dicho valor de abertura. Para dotar al modelo de una mayor similitud con el proceso real, se ha añadido una dinámica de primer orden con una constante de tiempo (𝜏𝑣) de unos pocos segundos: donde 𝐹𝑐,𝑟𝑒𝑞 es el caudal de refrigeración que se solicita con cada valor de abertura de la válvula. Nótese que al variar la abertura de la válvula (variable manipulable) se varía el caudal de refrigerante. De este modo, una nueva abertura (deseada) implicará un nuevo caudal (solicitado). Asimismo, la válvula modelada es de tipo “air-to- open”, es decir, la abertura de la misma aumenta de manera directamente proporcional a la señal de control. 𝜏𝑣𝑑𝐹𝑐 𝑑𝑡+𝐹𝑐=𝐹𝑐,𝑟𝑒𝑞 (12)
Modelado del reactor de tanque agitado continuo 6 Tabla 1. Datos de los parámetros del reactor CSTR con serpentín de refrigeración y condiciones de operación nominal Variable Condiciones de operación nominal Concentración de reactivo a la salida (𝒄𝑨) 55,50 𝑚𝑜𝑙/𝑚3 Temperatura del reactor (𝑻) 356,43 𝐾 (83,28 ℃) Temperatura de refrigeración a la salida (𝑻𝒄) 356,41 𝐾 (83,26 ℃) Caudal de reactivo/producto (𝑭) 7,50·10−4 𝑚3/𝑠 Caudal de refrigerante (𝑭𝒄) 1,50·10−3 𝑚3/𝑠 Concentración de reactivo a la entrada (𝒄𝑨𝒆) 930,00 𝑚𝑜𝑙/𝑚3 Temperatura a la entrada del reactor (𝑻𝒆) 315,00 𝐾 (41,85 ℃) Temperatura de refrigeración a la entrada (𝑻𝒄𝒆) 310,00 𝐾 (36,85 ℃) Volumen del reactor (V) 0,10 𝑚3 Término de transferencia de calor (𝑼𝑨) 4,88·104 𝑊/𝐾 Factor de frecuencia de la reacción (𝒌𝟎) 8,75·108 𝑠−1 Energía de activación (𝑬) 67,34·103 𝐽/𝑚𝑜𝑙 Calor de reacción (−𝚫𝑯) 8,37 𝑥 105 𝐽/𝑚𝑜𝑙 Densidad (𝝆) 1,00·103 𝑘𝑔/𝑚3 Calor específico (𝑪𝒑) 4,19·103 𝐽/(𝑘𝑔 𝐾) Peso específico del refrigerante (𝜸) 1,00 2.2 Modelo del reactor CSTR encamisado Suponiendo que la extracción de calor de la reacción exotérmica se lleva a cabo mediante una camisa de refrigeración (Figura 2) en lugar de a través de un serpentín, se asume lo siguiente: La pureza del reactivo a la entrada del reactor es del 100% La mezcla dentro del reactor está perfectamente agitada. El fluido de refrigeración en el interior de la camisa tiene una temperatura uniforme y la temperatura del caudal de refrigeración a la salida es igual a la temperatura en el interior de ésta. La densidad del reactivo a la entrada es igual a la densidad del producto a la salida. Existe un control de nivel en el tanque asociado al caudal de entrada, de forma que el volumen de la mezcla en el interior se mantiene constante en el tiempo. Se desprecian las pérdidas de calor hacia el exterior. Los calores específicos a presión constante (Cp) no cambian con la temperatura dentro del intervalo de operación, y coinciden en las corrientes del proceso a la entrada y salida del reactor.
7 Control de un reactor continuo de tanque agitado El balance de masa en el interior del reactor no varía con respecto al caso del reactor refrigerado mediante serpentín, por lo que se tiene: 𝑉𝑑𝑐𝐴 𝑑𝑡=𝐹(𝑐𝐴𝑒−𝑐𝐴)−𝑉𝑘𝑐𝐴 (13) El calor que se extrae del reactor por la camisa de refrigeración puede modelarse mediante la siguiente expresión: 𝑄=𝑈𝐴(𝑇−𝑇𝑐) (14) Por lo que la evolución de la temperatura en el reactor viene dada por: 𝑉𝜌𝐶𝑝𝑑𝑇 𝑑𝑡=𝜌𝐶𝑝𝐹(𝑇𝑒−𝑇)+(−Δ𝐻𝑅)𝑉𝑘𝑐𝐴+𝑈𝐴(𝑇𝑐−𝑇) (15) Y la evolución de la temperatura en la camisa de refrigeración puede expresarse como: 𝑑𝑇𝑐 𝑑𝑡=𝐹𝑐 𝑉𝑐(𝑇𝑐𝑒−𝑇𝑐)+𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝑉𝑐(𝑇−𝑇𝑐) (16) Las variables de entrada y salida del proceso son las mismas que en el caso del reactor con serpentín de refrigeración. Los parámetros de la camisa y las condiciones iniciales del modelo implementado en GUI se muestran en la Tabla 2. Nótese que el punto de equilibrio difiere respecto al caso del reactor refrigerado mediante serpentín (Tabla 1). Al igual que en el reactor con serpentín de refrigeración, existe una válvula de control localizada a la entrada de la camisa y modelada mediante la ecuación (12). Figura 2. Esquema de un reactor de tanque agitado continuo con camisa de refrigeración
Modelado del reactor de tanque agitado continuo 8 Tabla 2. Datos de los parámetros de la camisa de refrigeración y condiciones de operación nominal Variable Condiciones de operación nominal Concentración de reactivo a la salida (𝒄𝑨) 42,30 𝑚𝑜𝑙/𝑚3 Temperatura del reactor (𝑻) 360,96 𝐾 (87,81 ℃) Temperatura de refrigeración a la salida (𝑻𝒄) 355,16 𝐾 (82,01 ℃) Caudal de refrigerante (𝐅𝐜) 1,50·10−3𝑚3/𝑠 Temperatura de refrigeración a la entrada (𝑻𝒄𝒆) 310,00 𝐾 (36,85 ℃) Término de transferencia de calor (𝐔𝐀) 4,88· 104 𝑊/𝐾 Volumen de la camisa del reactor (𝐕𝐜) 0,01 𝑚3 Densidad ( 𝛒𝐜) 1,00· 103 𝑘𝑔/𝑚3 Calor específico ( 𝐂𝐜) 4,19·103 𝐽/(𝑘𝑔 𝐾) 2.3 Modelo del reactor CSTR encamisado con recirculación de refrigerante La adición de un lazo de recirculación de caudal en la camisa (Figura 3) permite obtener un mejor comportamiento, ya que minimiza los efectos de: Generación de puntos fríos y calientes, que conducen a una pobre homogeneización de la temperatura del fluido en la camisa de refrigeración. Tiempo muerto y ganancia del proceso que serían variables. Bajo coeficiente global de transferencia de calor a caudales bajos. Ensuciamiento a caudales bajos. Para el modelado de este proceso se tienen en cuenta las mismas consideraciones tomadas para el caso del reactor encamisado, además de las siguientes: Existe una bomba que mantiene constante el caudal de refrigeración en la camisa. En la zona de impulsión de la bomba, una vez derivado el caudal de recirculación, existe, en la purga de refrigerante, una válvula reguladora de presión aguas arriba de la misma. En la tubería de entrada a la camisa, aguas abajo del nodo de recirculación, existe una válvula que puede manipularse para variar el caudal de refrigerante fresco y así controlar la temperatura dentro del reactor.
9 Control de un reactor continuo de tanque agitado Figura 3. Esquema de un reactor de tanque agitado continuo encamisado con recirculación de refrigerante La evolución de la temperatura del caudal de refrigerante viene dada por la siguiente ecuación diferencial: 𝑑𝑇𝑐 𝑑𝑡=𝐹𝑐 𝑉𝑐(𝑇𝑗−𝑇𝑐)+𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝑉𝑐(𝑇−𝑇𝑐) (17) donde 𝑇𝑗 es la temperatura resultante de la mezcla entre el caudal de refrigeración de entrada al proceso y la fracción de flujo recirculado, procedente de la salida de la camisa. Realizando un balance de masa en los nodos se tiene: 𝐹𝑐=𝐹𝑐𝑒+𝐹𝑟𝑐 (18) 𝐹𝑐=𝐹𝑐𝑠+𝐹𝑟𝑐 (19) donde 𝐹𝑟𝑐 es el caudal recirculado y 𝐹𝑐𝑠 el caudal de purga a través de la válvula reguladora de presión. El balance de energía del proceso proporciona la siguiente expresión: 𝜌𝑐𝐹𝑐𝑒𝐶𝑝𝑐(𝑇𝑐𝑒−𝑇𝑟𝑒𝑓)+𝜌𝑐𝐹𝑟𝑐𝐶𝑝𝑐(𝑇𝑐−𝑇𝑟𝑒𝑓)=𝜌𝑐𝐹𝑐𝐶𝑝𝑐(𝑇𝑗−𝑇𝑟𝑒𝑓) (20) Tomando como temperatura de referencia 𝑇𝑟𝑒𝑓=273 𝐾, se tiene: 𝐹𝑐𝑇𝑗=𝐹𝑐𝑒𝑇𝑐𝑒+𝐹𝑟𝑐𝑇𝑐 (21) 𝑇𝑗=𝐹𝑐𝑒𝑇𝑐𝑒+𝐹𝑐𝑟𝑇𝑐 𝐹𝑐=𝑇𝑐+𝐹𝑐𝑒 𝐹𝑐(𝑇𝑐𝑒−𝑇𝑐) (22) Finalmente, sustituyendo 𝑇𝑗 en la ecuación (17), se tiene: 𝑑𝑇𝑐 𝑑𝑡=𝐹𝑐𝑒 𝑉𝑐(𝑇𝑐𝑒−𝑇𝑐)+𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝑉𝑐(𝑇−𝑇𝑐) (23)
Modelado del reactor de tanque agitado continuo 10 Nótese que este sistema de ecuaciones coincide con el descrito en el apartado 2.2, aunque esta nueva configuración da como resultado una mayor similitud entre el proceso real y el modelo matemático. Esto es debido, en parte, a una menor variación del coeficiente global de transferencia de calor, obtenido al mantener un caudal de refrigeración constante a través de la camisa. En la GUI, se han implementado las ecuaciones correspondientes al caso con recirculación, sin embargo, se ha optado por representar en la misma el esquema correspondiente a la Figura 2, con el objetivo de simplificar la información visual de la ventana principal. 2.4 Aproximación de retardos presentes en el proceso real Las ecuaciones anteriormente descritas han sido obtenidas a partir de simplificaciones y suposiciones iniciales, con el objetivo de obtener un equilibrio entre semejanza con la realidad y sencillez en la integración con métodos numéricos. Como consecuencia, las simulaciones realizadas no muestran el efecto de tiempo muerto característico de este tipo de procesos. Para corregir esta desviación en el modelado, es posible incluir en el sistema de ecuaciones diferenciales la dinámica de la pared del reactor, de manera que se proporciona una nueva capacidad térmica al modelo, que conduciría a un sistema de tercer orden. Sin embargo, para no aumentar con ello la complejidad del sistema, se ha optado por incluir un retraso en los datos de entrada del simulador. Este retraso puro se ha dimensionado para incluir también la no linealidad asociada a una mezcla imperfecta y al retardo debido a una posición alejada del sensor-transductor, así como al efecto debido a algunas de las hipótesis simplificadoras consideradas en el modelado. Además, para dotar al sistema de un mayor realismo, se ha dado un retraso diferente a cada variable interna del proceso (temperatura del reactor, concentración a la salida y temperatura de refrigeración), de forma que se obtienen tiempo muertos distintos para cada relación estado interno/entrada. 2.5 Linealización del modelo en torno al punto de funcionamiento lineal Para realizar un análisis de la controlabilidad del proceso se ha realizado una previa linealización del modelo del reactor en torno al punto de operación nominal. Para ello, se han calculado las matrices Jacobianas de la función vectorial multivariable que describe la dinámica del sistema con respecto al estado del mismo y a la entrada, resultando en una expresión en el espacio de estados. A continuación, se muestra el proceso de linealización del modelo de reactor CSTR encamisado. Partiendo del siguiente sistema de ecuaciones diferenciales: 𝑑𝑐𝐴 𝑑𝑡=𝑓1(𝑐𝐴,𝑇,𝑇𝑐,𝐹𝑐) (24) 𝑑𝑇 𝑑𝑡=𝑓2(𝑐𝐴,𝑇,𝑇𝑐,𝐹𝑐) (25) 𝑑𝑇𝑐 𝑑𝑡=𝑓3(𝑐𝐴,𝑇,𝑇𝑐,𝐹𝑐) (26) Se obtienen las matrices que describen la evolución del estado de un sistema lineal aproximado, a partir de las matrices Jacobianas:
11 Control de un reactor continuo de tanque agitado 𝐴= [ 𝜕𝑓1 𝜕𝑐𝐴𝜕𝑓1 𝜕𝑇 𝜕𝑓1 𝜕𝑇𝑐 𝜕𝑓2 𝜕𝑐𝐴𝜕𝑓2 𝜕𝑇 𝜕𝑓2 𝜕𝑇𝑐 𝜕𝑓3 𝜕𝑐𝐴𝜕𝑓3 𝜕𝑇 𝜕𝑓3 𝜕𝑇𝑐 ] | | 𝑐𝐴=𝑐𝐴0 𝑇=𝑇0 𝑇𝑐=𝑇𝑐0 𝐹𝑐=𝐹𝑐0 (27) 𝐵= [ 𝜕𝑓1 𝜕𝐹𝑐 𝜕𝑓2 𝜕𝐹𝑐 𝜕𝑓3 𝜕𝐹𝑐 ] | | 𝑐𝐴=𝑐𝐴0 𝑇=𝑇0 𝑇𝑐=𝑇𝑐0 𝐹𝑐=𝐹𝑐0 (28) Seguidamente, se muestran los sistemas obtenidos para cada modelo de reactor: Reactor CSTR con serpentín de refrigeración 𝐴= [ −𝐹𝑉−𝑘0𝑒−𝐸 𝑅𝑇0−𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇0 (−Δ𝐻𝑟) 𝜌𝐶𝑝𝑘0𝑒−𝐸 𝑅𝑇0−𝐹𝑉+(−Δ𝐻𝑟) 𝜌𝐶𝑝𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇0−𝜌𝑐𝐶𝑝𝑐𝐹𝑐0 𝑉𝐶𝑝𝜌(1−𝑒 −𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐0) ] (29) 𝐵= [ 0 𝜌𝑐𝐶𝑝𝑐 𝜌𝐶𝑝𝑉(𝑇0−𝑇𝑐𝑒)(𝑒−𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐0(1+ 𝑈𝐴 𝜌𝑐𝐶𝑝𝑐𝐹𝑐0)−1) ] (30) [Δ𝑐𝐴 Δ𝑇]=𝐴·[Δ𝑐𝐴 Δ𝑇]+𝐵·Δ𝐹𝑐 (31) Reactor CSTR encamisado 𝐴= [ −𝐹𝑉−𝑘0𝑒−𝐸 𝑅𝑇0−𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇00 (−Δ𝐻𝑟) 𝜌𝐶𝑝𝑘0𝑒−𝐸 𝑅𝑇0−𝐹𝑉−𝑈𝐴 𝑉𝜌𝐶𝑝+(−Δ𝐻𝑟) 𝜌𝐶𝑝𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇0𝑈𝐴 𝑉𝜌𝐶𝑝 0𝑈𝐴 𝑉𝑐𝜌𝑐𝐶𝑝𝑐 −𝐹𝑐0 𝑉𝑐−𝑈𝐴 𝑉𝑐𝜌𝑐𝐶𝑝𝑐 ] (32) 𝐵=[00 𝑇𝑐𝑒−𝑇𝑐0 𝑉𝑐] (33)
Modelado del reactor de tanque agitado continuo 12 [Δ𝑐𝐴 Δ𝑇 Δ𝑇𝑐]=𝐴·[Δ𝑐𝐴 Δ𝑇 Δ𝑇𝑐]+𝐵·Δ𝐹𝑐 (34) Reactor CSTR encamisado con recirculación del caudal de refrigeración: 𝐴= [ −𝐹𝑉−𝑘0𝑒−𝐸 𝑅𝑇0−𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇00 (−Δ𝐻𝑟) 𝜌𝐶𝑝𝑘0𝑒−𝐸 𝑅𝑇0−𝐹𝑉−𝑈𝐴 𝑉𝜌𝐶𝑝+(−Δ𝐻𝑟) 𝜌𝐶𝑝𝐸 𝑅𝑇02𝑐𝐴0𝑘0𝑒−𝐸 𝑅𝑇0𝑈𝐴 𝑉𝜌𝐶𝑝 0𝑈𝐴 𝑉𝑐𝜌𝑐𝐶𝑝𝑐 −𝐹𝑐𝑒0 𝑉𝑐−𝑈𝐴 𝑉𝑐𝜌𝑐𝐶𝑝𝑐 ] (35) 𝐵=[00 𝑇𝑐𝑒−𝑇𝑐0 𝑉𝑐] (36) [Δ𝑐𝐴 Δ𝑇 Δ𝑇𝑐]=𝐴·[Δ𝑐𝐴 Δ𝑇 Δ𝑇𝑐]+𝐵·Δ𝐹𝑐𝑒 (37) Una vez obtenidas las matrices y sustituidas las variables por los valores proporcionados en las Tablas 1 y 2, se ha analizado la controlabilidad del sistema lineal resultante, con el objetivo de comprobar si es posible el diseño de un controlador alrededor del punto de operación nominal del reactor. Dicho análisis se ha realizado comprobando el rango de la matriz de controlabilidad: 𝑀𝐶=[𝐵𝐴𝐵𝐴2𝐵… 𝐴𝑁1−1𝐵] (38) donde 𝑁1 es la dimensión del vector de estados, que tiene un orden de 𝑁1=2 en el caso del reactor con serpentín (temperatura del reactor y concentración de reactivo a la salida). En el reactor encamisado, se añade la dinámica de la temperatura de refrigeración, resultando en un vector de estados con un orden de 𝑁1=3. El resultado de la ecuación (38), para cada caso estudiado, coincide con la dimensión del correspondiente vector de estados 𝑁1, indicando que los sistemas son completamente controlables. Cabe destacar que el estudio de la controlabilidad de los modelos no lineales desarrollados en este capítulo es mucho más complejo, por lo que el resultado de este análisis tan sólo es válido de manera local y en el entorno cercano al punto nominal. Los límites de operación de las variables del proceso se recogen en la Tabla 3. Tabla 3. Rango de operación del reactor Variable Valor mínimo Valor máximo Concentración de reactivo a la salida (𝒄𝑨) 9,5 𝑚𝑜𝑙/𝑚3 74 𝑚𝑜𝑙/𝑚3 Temperatura del reactor (𝑻) 351,5 𝐾 (78,35 ℃) 387,8 𝐾 (114,65 ℃) Temperatura de refrigeración a la salida (𝑻𝒄) 345,6 𝐾 (72,45 ℃) 383,15 𝐾 (110 ℃) Caudal de refrigerante (𝐅𝐜) 7·10−4 𝑚3/𝑠 2·10−3 𝑚3/𝑠
13 3 SIMULACIÓN DEL PROCESO 3.1 Simulación de las ecuaciones mediante integración numérica Para integrar las ecuaciones del modelo en la interfaz gráfica se han implementado los métodos de integración numérica de Euler y Runge-Kutta de orden 4, debido que se tratan de métodos de paso fijo, de sencilla implementación y que, generalmente, proporcionan resultados precisos con modelos dinámicos simples. El algoritmo de Runge-Kutta es, generalmente, más preciso que el de Euler, aunque requiere un mayor esfuerzo computacional para su ejecución. Es por ello que se deja a elección del usuario el método utilizado por la aplicación para integrar las ecuaciones del modelo, a partir de una opción localizada dentro de la interfaz. A continuación, se detallan ambos algoritmos. Euler: Dada una ecuación diferencial expresada de la forma: 𝑑𝑦 𝑑𝑡=𝑓(𝑦,𝐻(𝑦,𝑢(𝑡)),𝑢(𝑡)) (39) Y dadas unas condiciones iniciales, puede obtenerse el valor en el instante siguiente: 𝑦𝑘+1=𝑦𝑘+Δ𝑡𝑘·𝑓(𝑦𝑘,𝐻(𝑦𝑘,𝑢𝑘),𝑢𝑘) (40) Runge-Kutta orden 4 Mediante este método se tiene: 𝑦1=𝑦𝑘 (41) 𝑢1=𝑢(𝑡𝑘) (42) 𝑦2=𝑦𝑘+12Δ𝑡𝑘𝑤1 (43) 𝑦3=𝑦𝑘+12Δ𝑡𝑘𝑤2 (44) 𝑦4=𝑦𝑘+Δ𝑡𝑘𝑤3 (45) 𝑢2=𝑢3=𝑢(𝑡𝑘+12Δ𝑡𝑘) (46) 𝑢4=(𝑡𝑘+Δ𝑡𝑘) (47) 𝑧𝑖=𝐻(𝑦𝑖,𝑢𝑖) (48) 𝑤𝑖=𝑓(𝑦𝑖,𝑧𝑖,𝑢𝑖) (49) 𝑦𝑘+1=𝑦𝑘+Δ𝑡𝑘·(16𝑤1+13𝑤2+13𝑤3+16𝑤4) (50)
Simulación del proceso 20 Para el guardado de datos, se tienen tres entradas: Una consiste en un indicador que se activa en caso de que el usuario seleccione la opción correspondiente al guardado de los datos. Las otras dos se corresponden con el número máximo de datos que se pueden guardar y con la orden de detención del proceso de grabación. Las salidas corresponden a las variables del proceso que son guardadas en cada ciclo de ejecución de la función “Callback” y un indicador que se activa en caso de haber finalizado con éxito la operación y, mediante el cual, la aplicación reconoce si existen o no datos guardados. Estas salidas se comunican como argumento de entrada a la ventana gráfica destinada a la identificación del proceso, que se detalla en profundidad en el siguiente capítulo. Con respecto al almacenamiento del historial, los datos comienzan a registrarse desde el inicio de la simulación y, una vez alcanzado el número máximo de datos establecido, se elimina el primer valor registrado de cada variable y se almacena el siguiente, manteniendo constante, a partir de este momento, el tamaño de la base de datos.
21 4 IDENTIFICACIÓN DEL PROCESO ebido a que el sistema es no lineal, se ha incluido en la interfaz gráfica una herramienta que permite obtener modelos lineales aproximados del sistema, realizando previamente una simulación de ensayo con entrada en escalón. Estos modelos lineales son de gran utilidad para obtener un ajuste inicial del controlador en el punto de operación deseado, utilizando reglas de sintonización a partir de los parámetros identificados. Por otro lado, también resulta útil para el diseño de controladores más avanzados, como es el caso del control predictivo, en el que es necesario un modelo de predicción, o del control adaptativo “Gain Scheduling”, en el que resulta conveniente conocer distintos modelos lineales que se ajusten adecuadamente a cada punto de operación del proceso. En este capítulo, se sientan las bases de la identificación paramétrica de sistemas dinámicos en tiempo discreto a partir de datos obtenidos mediante ensayos experimentales. Luego, se realiza una descripción sobre la herramienta de identificación implementada en la interfaz gráfica, así como los métodos empleados para la detección y prevención de errores en el proceso de estimación. 4.1 Identificación paramétrica de sistemas lineales aproximados La identificación paramétrica abarca un amplio conjunto de métodos utilizados para la estimación de modelos matemáticos de procesos físicos a partir de datos obtenidos experimentalmente, ya sea de ensayos reales realizados en planta o a través de simulación. A diferencia de los métodos de identificación no paramétricos, este proceso da como resultado, de forma directa, los valores de los parámetros de la estructura matemática seleccionada para tal fin. Generalmente, estos parámetros se obtienen a partir de la definición y resolución de un problema de optimización sujeto a un determinado criterio, por lo que es este carácter óptimo lo que permite la obtención de modelos más fidedignos que los proporcionados por los clásicos métodos no paramétricos. El proceso general de identificación puede resumirse en los siguientes apartados: 1. Realización del ensayo en planta o ejecución de la simulación y registro de datos. 2. Preprocesamiento de los datos. 3. Selección del modelo matemático a identificar. 4. Selección del método de identificación. 5. Obtención de los parámetros óptimos, de acuerdo al criterio de optimización. 6. Validación del modelo. Este ensayo experimental suele consistir en la excitación de la planta a partir de cambios de tipo escalón en la entrada del sistema, de manera que, partiendo de un punto de operación nominal, se obtenga información que permita captar la dinámica del proceso. El muestreo de cada tipo de dato necesario se realiza de forma simultánea y sincronizada, permitiendo, posteriormente, relacionar en el tiempo un dato de entrada con el correspondiente valor de salida. El modelo obtenido mediante el proceso es normalmente lineal en tiempo discreto, de forma que se simplifica el diseño de controladores en el entorno cercano al punto de operación alrededor del cual se ha llevado a cabo el experimento. D
Identificación del proceso 22 Los tipos de modelos dinámicos más utilizado para la identificación paramétrica se detallan a continuación [6]: Modelo autorregresivo y exógeno (ARX): donde 𝜉(·) representa un ruido blanco con media cero y varianza 𝜎2. Se observa que la parte autorregresiva está formada por la secuencia de términos de salida en el pasado y la parte exógena por los términos de la secuencia de entrada. 𝑛𝑎 y 𝑛𝑏 son, respectivamente, el orden de la parte autorregresiva y de la parte exógena [6]. Nótese que las variables de entrada y salida son incrementales respecto a una posición de equilibrio. La función de transferencia resultante, en el dominio de 𝒵 es la siguiente: 𝑦(𝑧)=𝑏1𝑧−1+𝑏2𝑧−2+⋯+𝑏𝑛𝑏𝑧−𝑛𝑏 1+𝑎1𝑧−1+𝑎2𝑧−2+⋯+𝑎𝑛𝑎𝑧−𝑛𝑎u(z) + 1 1+𝑎1𝑧−1+𝑎2𝑧−2+⋯+𝑎𝑛𝑎𝑧−𝑛𝑎𝜉(𝑧) (54) 𝑦(𝑧)=𝐵(𝑧) 𝐴(𝑧)𝑢(𝑧)+1 𝐴(𝑧)𝜉(𝑧) (55) Puede comprobarse que el caso ARX equivale a una dinámica lineal a la que se le añade un ruido blanco filtrado por su propia ecuación característica (Figura 5). Figura 5. Diagrama de bloques de un modelo con estructura ARX Para la identificación, los parámetros de ajuste son: 𝜃=[𝑎1,𝑎2,…,𝑎𝑛𝑎,𝑏1,𝑏2,…,𝑏𝑛𝑏] (56) Modelo autorregresivo y exógeno de media móvil (ARMAX): 𝑦(𝑡)=𝑎1𝑦(𝑡−1)+𝑎2𝑦(𝑡−2)+⋯+𝑎𝑛𝑎𝑦(𝑡−𝑛𝑎)+ + 𝑏1𝑢(𝑡−1)+𝑏2𝑢(𝑡−2)+⋯+𝑏𝑛𝑏(𝑡−𝑛𝑏)+𝜐(𝑡) (57) Siendo 𝜐(𝑡) el término de media móvil asociado al ruido: 𝜐(𝑡)=𝜉(𝑡)+𝑐1𝜉(𝑡−1)+⋯+𝑐𝑛𝑐𝜉(𝑡−𝑛𝑐) (58) 𝑦(𝑡)=𝑎1𝑦(𝑡−1)+𝑎2𝑦(𝑡−2)+⋯+𝑎𝑛𝑎𝑦(𝑡−𝑛𝑎) + 𝑏1𝑢(𝑡−1)+𝑏2𝑢(𝑡−2)+⋯+𝑏𝑛𝑏𝑢(𝑡−𝑛𝑏)+𝜉(𝑡) (53)
23 Control de un reactor continuo de tanque agitado donde 𝑛𝑐 es el orden de la parte de media móvil [6]. Se observa que, a diferencia del caso ARX, se consideran los valores de la señal de ruido blanco en el pasado, por lo que se tiene un modelo con un mayor número de parámetros, dificultando el proceso de identificación, pero obteniendo, generalmente, mejores resultados. La función de transferencia resultante, en el dominio de 𝒵 es la siguiente: 𝑦(𝑧)=𝑏1𝑧−1+𝑏2𝑧−2+⋯+𝑏𝑛𝑏𝑧−𝑛𝑏 1+𝑎1𝑧−1+𝑎2𝑧−2+⋯+𝑎𝑛𝑎𝑧−𝑛𝑎u(z) +1+𝑐1𝑧−1+𝑐2𝑧−2+⋯+𝑐𝑛𝑐𝑧−𝑛𝑐 1+𝑎1𝑧−1+𝑎2𝑧−2+⋯+𝑎𝑛𝑎𝑧−𝑛𝑎𝜉(𝑧) (59) 𝑦(𝑧)=𝐵(𝑧) 𝐴(𝑧)𝑢(𝑧)+𝐶(𝑧) 𝐴(𝑧)𝜉(𝑧) (60) A continuación, se muestra el diagrama de bloques de la estructura ARMAX (Figura 6): Figura 6. Diagrama de bloques de un modelo con estructura ARMAX Para la identificación, los parámetros de ajuste son: 𝜃=[𝑎1,𝑎2,…,𝑎𝑛𝑎,𝑏1,𝑏2,…,𝑏𝑛𝑏,𝑐1,𝑐2,…,𝑐𝑛𝑐] (61) Modelo Box y Jenkins: Este tipo de modelo se trata de un caso general, que engloba a otros tipos de modelos, cuya expresión en el dominio de 𝒵 viene dado por: 𝑦(𝑧)=𝐵(𝑧) 𝐴(𝑧)𝑢(𝑧)+𝐶(𝑧) 𝐷(𝑧)𝜉(𝑧) (62)
Identificación del proceso 24 A continuación, se muestra de forma esquematizada un modelo Box y Jenkins (Figura 7): Figura 7. Diagrama de bloques de un modelo con estructura Box and Jenkins Para la identificación, los parámetros de ajuste son: 𝜃=[𝑎1,𝑎2,…,𝑎𝑛𝑎,𝑏1,𝑏2,…,𝑏𝑛𝑏,𝑐1,𝑐2,…,𝑐𝑛𝑐,𝑑1,𝑑2,…,𝑑𝑛𝑑] (63) Nótese que los órdenes de la parte autorregresiva y la parte exógena pueden seleccionarse para obtener un orden determinado en la función de transferencia de la dinámica identificada. En el control automático, suelen utilizarse para el diseño de controladores sistemas de primer y segundo orden: Sistemas de primer orden con tiempo muerto En tiempo continuo y en el dominio de ℒ, la expresión general de un sistema de primer orden con tiempo muerto viene dada por la siguiente función de transferencia: 𝑦(𝑠) 𝑢(𝑠)=𝐾𝑝 𝜏𝑝𝑠+1𝑒−𝜏𝑚𝑝𝑠 (64) donde 𝐾𝑝 es la ganancia estática, 𝜏𝑝 es la constante de tiempo del sistema y 𝜏𝑚𝑝 el tiempo muerto. Empleando un retenedor de orden cero (ZOH), la equivalencia en tiempo discreto puede expresarse mediante la siguiente relación de recurrencia: 𝑦(𝑘)=−𝑎·𝑦(𝑘−1)+𝑏·𝑢(𝑘−1−𝑛𝑑𝑒𝑙) (65) donde 𝑛𝑑𝑒𝑙 es el número de pasos que tarda un cambio en la entrada en afectar a la señal de salida. Sistemas de segundo orden con tiempo muerto La función de transferencia en tiempo continuo que describe la dinámica de un sistema de segundo orden con tiempo muerto puede expresarse mediante la siguiente ecuación: 𝑦(𝑠) 𝑢(𝑠)=𝐾𝑝𝜔𝑛2 𝑠2+2𝜁𝜔𝑛𝑠+𝑤𝑛2𝑒−𝜏𝑚𝑝𝑠 (66) donde 𝜔𝑛 es la frecuencia natural del proceso y 𝜁 la constante de amortiguamiento.
25 Control de un reactor continuo de tanque agitado Empleando un retenedor de orden cero (ZOH), la equivalencia en tiempo discreto puede expresarse mediante la siguiente relación de recurrencia: 𝑦(𝑘)=−𝑎1·𝑦(𝑘−1)−𝑎2·𝑦(𝑘−2)+𝑏1·𝑢(𝑘−𝑛𝑑𝑒𝑙−1) +𝑏2·𝑢(𝑘−𝑛𝑑𝑒𝑙−2) (67) Asimismo, para cualquier estructura matemática descrita en este apartado, los datos de la parte exógena pueden retrasarse 𝑛𝑑𝑒𝑙 unidades para incluir en la identificación el efecto de tiempo muerto. 4.2 Método de identificación paramétrica por mínimos cuadrados El método de identificación implementado en la interfaz gráfica es el de mínimos cuadrados, cuyo procedimiento se detalla a continuación. Dado un conjunto de datos de entrada y salida obtenidos por simulación o mediante ensayos reales, se define como error de modelado: 𝐸=𝑌−𝑋θ (68) donde 𝑋 es la matriz que contiene los datos de entrada, 𝑌 los de salida y 𝜃 es el vector paramétrico. El producto 𝑋𝜃 representa, por tanto, la estimación de la salida obtenida a partir de los datos de entrada y de los parámetros a estimar. Nótese que el conocimiento de los valores pasados de las entradas y salidas no proporcionan ninguna información útil para evaluar 𝜉(𝑡) y, por ello, esta no se incluye en la expresión del estimador [6]. Se debe minimizar la suma de los errores al cuadrado, por lo que puede definirse el siguiente problema de optimización: min 𝐸𝐸𝑇𝐸=min θ(𝑌−𝑋θ)𝑇(𝑌−𝑋𝜃) (69) Este problema tiene solución explícita dada por: 𝜃=(𝑋𝑇𝑋)−1𝑋𝑇𝑌 (70) Debido a que los reactores descritos en el apartado 2 se tratan de modelos de reducida complejidad, se ha optado por la selección de una estructura ARX para el algoritmo de identificación implementado en la GUI, con órdenes 𝑛𝑎=1 y 𝑛𝑏=1, que equivale a un sistema de primer orden (ecuación (65)). Nótese que un valor de 𝑛𝑎=1 implica que la salida en el instante anterior, 𝑦𝑘−1, tiene repercusión en la salida actual del proceso, 𝑦𝑘. Por otro lado, también puede comprobarse que la estructura ARX definida mediante la ecuación (53), incluye el retraso debido al tiempo de muestreo, al no considerarse el efecto de 𝑢𝑘. Con el objetivo de tener en cuenta el tiempo muerto en la estimación por mínimos cuadrados, se ha considerarlo dicho parámetro como una información conocida a priori, estableciendo un retraso de 𝑛𝑑𝑒𝑙 unidades de Δ𝑡𝑘 en los datos exógenos. El vector paramétrico, en este caso, presenta la siguiente estructura: 𝜃=[𝑎𝑏] (71)
Identificación del proceso 26 donde 𝑎 y 𝑏 son los parámetros del modelo lineal de primer orden (ecuación (65)). Por lo que debe seleccionarse una base de datos de entrada que contenga los valores de 𝑦(𝑘−1) y 𝑢(𝑘−1−𝑛𝑑𝑒𝑙) en cada instante de tiempo y una base de datos de salida con los valores de 𝑦(𝑘). Las filas de las matrices 𝑋 e 𝑌 son, por tanto: 𝑋𝑘,∗=[−𝑦𝑘−1 𝑢𝑘−1−𝑛𝑑𝑒𝑙] (72) 𝑌𝑘,∗=𝑦𝑘 (73) donde 𝑢 e 𝑦 son, respectivamente, los vectores de datos de entrada y de salida incrementales respecto al punto de equilibrio inicial. En la GUI, una vez obtenidos los datos experimentales, se aplica un filtro de media móvil sobre la base de datos, con el objetivo de reducir el impacto del ruido de medida sobre la identificación. Los valores de salida del filtro se utilizan como una nueva base de datos para el algoritmo de identificación paramétrica. En la Figura 8 puede observarse el suavizado obtenido tras aplicar el filtro a los datos medidos en un ensayo con entrada en escalón. Figura 8. Suavizado de datos mediante filtro de media móvil Una vez realizado el proceso de identificación desde la aplicación, se añade un valor aleatorio, dentro de un rango estrecho de valores, a cada parámetro obtenido, incorporando así el efecto potencial de las incertidumbres presentes en el experimento con el proceso real. Los números aleatorios son generados a partir de una distribución normal de media cero y con una desviación típica que depende del tipo de modelo seleccionado, del tipo de entrada (variable manipulada o perturbación) y del parámetro en cuestión. En la Figura 9, se muestra una comparación entre la respuesta del proceso modelado y la del sistema lineal identificado para las mismas condiciones de entrada a la simulación. El experimento se ha realizado variando el valor de abertura de la válvula de control un 5 % respecto a su posición de equilibrio cada 10 minutos. Puede comprobarse que existen pequeñas discrepancias entre las respuestas debido a las no linealidades del proceso.
27 Control de un reactor continuo de tanque agitado Figura 9. Comparación entre proceso modelado y proceso identificado empíricamente para distintos modelos de reactor CSTR En el Anexo B se detallan los datos de simulación, así como los resultados numéricos obtenidos en este experimento. 4.3 Transformación del modelo de tiempo discreto a tiempo continuo El método de identificación paramétrica descrito en el apartado anterior da como resultado un modelo lineal en tiempo discreto. Sin embargo, es importante destacar que la mayoría de métodos analíticos clásicos para el diseño de controladores, requieren de una expresión en tiempo continuo. Por ello, el algoritmo implementado realiza una transformación del sistema identificado para obtener una función de transferencia continua. Esta transformación puede llevarse a cabo mediante varios métodos. Los más utilizados son: Retenedor de orden cero (“Zero-Order Hold”): Se asume que la entrada es constante entre periodos de muestreo. Retenedor de primer orden (“First-Order Hold”): Se asume que la entrada es lineal creciente o decreciente entre periodos de muestreo. Transformación bilineal o método de Tustin: Aproximación de la derivada.
Identificación del proceso 28 El mantenedor de orden cero resulta muy útil para obtener una buena aproximación de la respuesta temporal para entradas en escalón, por lo que se ha optado por implementar este método en la interfaz gráfica. La transformación de tiempo continuo a discreto [7] viene dado por: 𝐺(𝑧)=(1−𝑧−1)𝒵{𝐺(𝑠) 𝑠} (74) donde para un sistema de primer orden, se tiene: 𝒵{𝑘𝑝 𝑠+𝑝 𝑠}=𝑘 (1−𝑒−𝑝Δ𝑡𝑘)𝑧 (𝑧−1)(𝑧−𝑒−𝑝Δ𝑡𝑘) (75) 𝐺(𝑧)=𝑘1−𝑒−𝑝Δ𝑡𝑘 𝑧−𝑒−𝑝Δ𝑡𝑘=𝑘(1−𝑒−𝑝Δ𝑡𝑘)𝑧−1 1−𝑒−𝑝Δ𝑡𝑘·𝑧−1 (76) Los modelos identificados por la aplicación vienen dados por la siguiente expresión: 𝐺(𝑧)=𝑏𝑧−1 1+𝑎𝑧−1 (77) Con lo que puede obtenerse la correspondencia con la constante de tiempo de la siguiente forma: 𝜏=1𝑝=Δ𝑡𝑘 |−ln(−𝑎)| (78) donde el parámetro 𝑎 se obtiene de la identificación paramétrica, y tiene valor negativo si se considera una base de datos acorde a la estructura representada en las ecuaciones (72) y (73). La ganancia estática puede obtenerse aplicando el teorema del valor final para sistemas en tiempo discreto: 𝑦(∞)=lim 𝑘→∞[𝑦(𝑘)]=lim 𝑧→1[(1−𝑧−1)𝑦(𝑧)] (79) Aplicando un escalón unitario se tiene: 𝑦(∞)=lim 𝑧→1[𝑏𝑧−1 1+𝑎𝑧−1]=𝑏 1+𝑎 (80) donde 𝑏 es el otro elemento del vector paramétrico obtenido en la identificación. El tiempo muerto debe hacerse igual al tiempo de retraso: 𝜏𝑚𝑝=𝑛𝑑𝑒𝑙·Δ𝑡𝑘 (81) Nótese que el valor de 𝑛𝑑𝑒𝑙 se considera como información conocida y pertenece al conjunto de argumentos de entrada del algoritmo. Esto implica que el tiempo muerto y la constante de tiempo no guardan ninguna relación desde el punto de vista de la identificación. Finalmente, se obtiene la siguiente expresión: 𝐺(𝑠)=𝑏 1+𝑎 Δ𝑡𝑘 |−ln(−𝑎)|𝑠+1𝑒−𝜏𝑚𝑝𝑠 (82)
29 Control de un reactor continuo de tanque agitado 4.4 Herramienta de ensayos con entrada en escalón para identificación de modelos aproximados La interfaz gráfica dispone de una ventana dedicada a la identificación de modelos lineales a partir de ensayos con entrada en escalón, de forma que el usuario puede introducir los valores de abertura de la válvula (variable manipulada) y los instantes en los que éstos se van a aplicar para, posteriormente, simular el proceso y realizar una identificación paramétrica por mínimos cuadrados con los resultados obtenidos. En la Figura 10 se ilustra con un diagrama de bloques el sistema lineal completo: Figura 10. Diagrama de bloques del sistema lineal completo En el proceso real, el ensayo con entrada en escalón consiste en cambiar la variable manipulada estando el controlador en manual. Este dispositivo dispondría de un botón con una escala porcentual, correspondiente al % de la abertura deseada para la válvula de control. Esta señal porcentual será convertida internamente a una señal eléctrica digital en tensión, y llegará al actuador eléctrico del elemento final de control en forma de señal eléctrica continua (𝑚𝐴,𝑜 𝑉) mediante una convertidor digital-analógico asociado al propio controlador. El actuador de la válvula convertirá la señal de 𝑚𝐴 (𝑜 𝑉) a %, y luego, esa señal porcentual vinculada a la posición del vástago de la válvula hará que cambie la posición de su obturador (es el mismo cambio por ser solidarios vástago y obturador), lo cual cambiará el caudal de fluido refrigerante, siguiendo una dinámica en la que la ganancia estática será la sensibilidad de la válvula (𝑙/𝑚𝑖𝑛 / %). A continuación, una vez cambia el caudal, cambiará el valor de la variable de interés en el proceso propiamente dicho (℃), con una dinámica que será la más lenta de todas las dinámicas implicadas en el ensayo. El valor de la variable controlada solo puede conocerse midiéndola, es decir, mediante un sensor-transmisor, cuya función de transferencia es la sensibilidad de este (𝑚𝐴 𝑜 𝑉/℃), pues se ha asumido que su dinámica es despreciable. La señal de 𝑚𝐴 (𝑜 𝑉) llegará al controlador, que la transformará en una señal digital de tensión mediante un convertidor analógico-digital (𝑉/(𝑚𝐴 𝑜 𝑉)), que luego podría visualizarse en el display del propio controlador como una señal porcentual (%/𝑉) o ingenieril (℃/𝑉). La Figura 10 recoge dos bloques asociados a la válvula, que son 𝐺𝑎𝑣 y 𝐺𝑣, correspondientes al actuador y al cuerpo de la válvula de control. El sensor-transmisor es necesario, ya que es la forma de saber cuál es el valor de la variable controlada en un proceso real que tiene lugar en una planta, y forma parte de la identificación. La cadena completa válvula de control + proceso propiamente dicho + sensor-transmisor conforma la dinámica global con la que deberá sintonizarse el controlador, cuya ganancia será adimensional. En la identificación implementada en la GUI se registran los datos del ensayo mediante el sensor-transmisor que transmite la señal en una magnitud eléctrica, posteriormente transformada a ℃ para la salida del proceso. Por otro lado, se registran los datos de abertura de válvula introducidos por el usuario que, en el proceso real, se obtendrían mediante el propio controlador configurado en manual, que puede tener una salida a un dispositivo de adquisición de datos, el cual también registraría la señal de temperatura. Con estos valores se lleva a cabo el proceso de identificación, dando lugar a una función de transferencia con unidades de ganancia de [℃/%], lo que permite desacoplar la identificación de las variables eléctricas. Tomando como base el diagrama de la Figura 10, se incorporan dos bloques conversores, 𝑓𝑠𝑡 y 𝑓𝑐, de modo que hay conversiones adicionales que afectan a las señales eléctricas: a ℃ para la salida del proceso y de % para la abertura de la válvula. Esto implica usar convertidores analógico-digital + digital a magnitud física (en ℃) (𝑓𝑠𝑡 global dado en ℃/(𝑉 𝑜 𝑚𝐴)) y porcentual a digital + digital-analógico (𝑓𝑐 global dado en (𝑉 𝑜 𝑚𝐴)/%), respectivamente. Por tanto, el producto de las funciones de transferencia es el siguiente:
Sistema de control por realimentación 36 Nótese que esta última expresión es válida para la regulación del proceso, pero no para su control sobre distintos puntos de operación, ya que no admite ningún valor de consigna para el seguimiento de referencias. Acción integral: Corresponde al término de la señal de control que es proporcional a la integral del error. Un controlador proporcional-integral (PI) para una estructura estándar (algoritmo ideal) tiene la siguiente expresión: 𝑢(𝑡)=𝐾𝑐(𝑒𝑟(𝑡)+1 𝑇𝑖∫𝑒𝑟(𝜏)𝑑𝜏 𝑡 0) (86) Mediante la integración, se proporciona información al controlador sobre el error acumulado a lo largo del tiempo, de forma que, si el sistema se estabiliza, el error de seguimiento tiende a ser nulo para valores finitos del tiempo integral (error en régimen permanente nulo). Este parámetro también afecta a la respuesta transitoria, obteniéndose una evolución más agresiva para valores bajos de 𝑇𝑖 o una mayor robustez para valores altos. Como contrapartida, la acción integral introduce un polo en el origen, lo que conlleva a una pérdida de estabilidad relativa. Esto puede reflejarse en una respuesta oscilatoria que, generalmente, se considera indeseable. Cuando un controlador incorpora acción integral, existe la posibilidad de que ocurra el denominado efecto “windup”. Este fenómeno tiene lugar cuando los actuadores saturan debido a sus limitaciones físicas, provocando una continua acumulación del error. Si esta saturación se alcanza debido a un cambio significativo en el valor nominal de consigna del controlador o de la perturbación, es posible que la integral del error aumente o disminuya indefinidamente, de modo que deje al sistema fuera de control incluso más allá del momento en el que la excesiva perturbación o referencia cambien, tardando mucho más tiempo en conseguir mantener bajo control el sistema, dando lugar a oscilaciones o incluso pudiendo llegar a inestabilizar al sistema. Para evitarlo, existen algunos métodos “anti-windup” en tiempo continuo, desarrollados para controladores analógicos, como es el caso de algoritmos incrementales o “back calculation”, sin embargo, debido a que, actualmente, la mayoría de controladores son digitales, suelen emplearse algoritmos en tiempo discreto. Una solución simple consiste en desactivar la acción integral cuando esta se satura, permitiendo su activación una vez que la saturación del (o de los) elemento(s) de control desaparece. Acción derivativa: Corresponde al término de la señal de control que es proporcional a la derivada del error. Un controlador proporcional-derivativo (PD) para una estructura estándar tiene la siguiente expresión: 𝑢(𝑡)=𝐾𝑐(𝑒𝑟(𝑡)+𝑇𝑑𝑑𝑒𝑟(𝑡) 𝑑𝑡 ) (87) La acción derivativa proporciona al controlador información sobre el valor de la pendiente de la curva del error, de forma que se obtiene una predicción de su evolución, permitiendo mejorar la estabilidad y reducir la sobreoscilación de la respuesta temporal. Nótese que la acción derivativa pura no es físicamente realizable, al no poder obtenerse la derivada del error en el instante actual (tangente de la curva del error en el último punto calculado), por lo que deben llevarse a cabo aproximaciones de la misma.
37 Control de un reactor continuo de tanque agitado La principal desventaja de derivar la señal de error es que puede producirse un “pico derivativo” (“derivative kick”). Este fenómeno ocurre cuando se realiza un cambio brusco en el valor de consigna, produciéndose una rápida variación en la derivada del error. El resultado es un repentino aumento del esfuerzo de control que puede producir efectos indeseados en la salida de la planta. Además, debido a que el ruido de medida induce variaciones abruptas sobre la señal de error, se producen cambios indeseados en su derivada que son amplificados por la acción derivativa. Por ello, no debe usarse esta acción en los procesos donde el ruido es alto (dando lugar a un controlador PI). Estas desventajas se mitigan parcialmente empleando un filtro paso bajo, que viene dado por la siguiente ecuación. 𝛼𝑇𝑑𝑑𝑢𝐷(𝑡) 𝑑𝑡 +𝑢𝐷(𝑡)=𝐾𝑐𝑇𝑑𝑑𝑒𝑟(𝑡) 𝑑𝑡 (88) donde 𝑢𝐷(𝑡) es la contribución de la acción derivativa a la señal de control. El parámetro 𝛼 se relaciona con la constante de tiempo del filtro mediante la siguiente expresión: 𝜏𝑓=𝛼𝑇𝑑 (89) Realizando la transformada de Laplace a la ecuación (88) se tiene: 𝑢𝐷(𝑠) 𝑒𝑟(𝑠)=𝐾𝑐𝑇𝑑𝑠 𝛼𝑇𝑑𝑠+1 (90) Para reducir el pico derivativo, la acción derivativa debe actuar sobre la variable controlada, resultando en la siguiente expresión para el controlador PD: 𝑢(𝑡)=𝐾𝑐(𝑒𝑟(𝑡)−𝑇𝑑𝑑𝑦(𝑡) 𝑑𝑡 ) (91) La estructura de control PID implementada en la interfaz gráfica consiste en el algoritmo estándar. La expresión del controlador resultante se indica a continuación: 𝑢(𝑡)=𝐾𝑐(𝑒𝑟(𝑡)+1 𝑇𝑖∫𝑒𝑟(𝜏)𝑑𝜏 𝑡 0)+𝑢𝐷(𝑡) (92) donde 𝑢𝐷(𝑡) depende de la configuración aportada por el usuario. Expresando la ecuación en el dominio de ℒ para el caso de actuación sobre el error, se tiene: 𝑢(𝑠)=𝐾𝑐(1+1 𝑇𝑖𝑠+𝑇𝑑𝑠 𝛼𝑇𝑑𝑠+1)𝑒𝑟(𝑠) (93) El usuario, desde una ventana de configuración, puede: Incluir o eliminar la acción integral y/o la acción derivativa del algoritmo de control PID. Sintonizar un control proporcional puro con actuación sobre la señal controlada o señal de error. Configurar o eliminar el filtro de la acción derivativa. Seleccionar si la acción derivativa actúa sobre la variable controlada o sobre la señal de error.
Sistema de control por realimentación 38 5.1.2 Discretización del controlador PID Para poder llevar a cabo una simulación por computación del proceso controlado, es necesario discretizar la ecuación del controlador PID, de forma que el valor de la señal de control se actualice en cada ciclo del bucle descrito en la Tabla 4 (Algoritmo 1: Simulación del sistema controlado). Debido a su reducido coste computacional y a su adecuada precisión, se ha utilizado el método de Euler para aproximar la integración de la señal de error: 𝐼𝐸𝑘+1=𝐼𝐸𝑘+𝐸𝑘·Δ𝑡𝑘 (94) La derivada de error puede aproximarse mediante el método de Euler hacia atrás: 𝐷𝑒𝑟,𝑘=𝑒𝑟,𝑘−𝑒𝑟,𝑘−1 Δ𝑡𝑘 (95) De esta forma, la salida del controlador sin filtro en la acción derivativa en cada instante viene dada por: 𝑢𝑘=𝑢0+𝐾𝑐𝑑(𝑒𝑟,𝑘+1 𝑇𝑖𝑑·𝐼𝐸𝑘+𝑇𝑑𝑑·𝐷𝐸𝑘) (96) Para aproximar el término derivativo, es preferible analizar las ecuaciones en el dominio de 𝒵. La transformada 𝒵 de una secuencia temporal viene dado por: 𝒵{𝑥(𝑛)}=∑−𝑥(𝑘)𝑧−𝑘 ∞ 𝑘=−∞ (97) donde la relación entre las variables complejas 𝑧 y 𝑠, se muestra a continuación: 𝑧=𝑒Δ𝑡𝑘·𝑠 (98) Nótese que el inverso de 𝑧 equivale a un retardo de Δ𝑡𝑘 unidades de tiempo, por lo que, a partir de la ecuación (95), se deduce la expresión en el dominio de 𝒵 de la derivada de la señal de salida: 𝐷𝐸(𝑧)=1−𝑧−1 Δ𝑡𝑘𝑒𝑟(𝑧) (99) Existen diversos métodos para aproximar el filtrado de la acción derivativa en tiempo discreto. Debido a su simple implementación y a su reducida amplificación de ruido respecto a otros métodos, se ha optado por sustituir el operador de Laplace “s” por una diferencia regresiva [10]: 𝑠=1−𝑧−1 Δ𝑡𝑘 (100) Aplicando la ecuación (100) en el término derivativo del controlador PID en tiempo continuo, dado por la ecuación (90) se tiene: 𝑢𝐷(𝑧) 𝑒𝑟(𝑧)=𝐾𝑐𝑇𝑑𝑑𝑁𝑓−𝑇𝑑𝑑𝑁𝑓𝑧−1 𝑁𝑓Δ𝑡𝑘+𝑇𝑑𝑑−𝑇𝑑𝑑𝑧−1 (101)
39 Control de un reactor continuo de tanque agitado donde 𝑇𝑑𝑑 es el tiempo derivativo y 𝑁𝑓 es un parámetro de diseño del filtro, que se corresponde con el inverso de 𝛼, y cuya relación con la constante de tiempo del mismo viene dada por: 𝜏𝑓=𝑇𝑑𝑑 𝑁𝑓 (102) Expresando la ecuación (101) mediante una relación de recurrencia se tiene: 𝑢𝐷𝑘=𝑇𝑑𝑑 𝑇𝑑𝑑+𝑁𝑓Δ𝑡𝑘𝑢𝐷𝑘−1+𝐾𝑐𝑇𝑑𝑑𝑁𝑓 𝑇𝑑𝑑+𝑁𝑓Δ𝑡𝑘(𝑒𝑟,𝑘−𝑒𝑟,𝑘−1) (103) La Figura 13 muestra el desempeño de la acción derivativa para tres casos con diferente de implementación. Los resultados se han obtenido realizando varias simulaciones del reactor CSTR encamisado para las mismas condiciones del ensayo, realizando un cambio en la referencia en el instante 𝑡=0,5 𝑚𝑖𝑛, con un tiempo derivativo de 𝑇𝑑𝑑=0,1 y sin simular ruido de medida. Puede comprobarse que derivando la señal de salida en lugar del error y añadiendo un filtro se obtiene una mayor suavidad en la señal de control y se evita el efecto del pico derivativo. Figura 13. Comparación de desempeño entre distintos tipos de acción derivativa Con la adición de este filtro, la salida del controlador en cada instante se obtiene con la siguiente expresión: 𝑢𝑘=𝑢0+𝐾𝑐𝑑(𝐸𝑘+1 𝑇𝑖𝑑·𝐼𝐸𝑘)+𝑢𝐷𝑘 (104)
Sistema de control por realimentación 40 donde para pasar de variables incrementales a variables absolutas, se ha añadido el término 𝑢0, que corresponde al valor nominal de la abertura de la válvula de control. Para evitar el efecto “windup”, se ha optado por desactivar la acción integral en caso de saturación de la válvula de control dispuesta en la línea de refrigerante. 5.1.3 Implementación en la GUI del algoritmo de control PID Debido a que el controlador simulado en la GUI es digital, se han incluido módulos de conversión, con el objetivo de procesar los datos de los sensores y obtener una ganancia estática adimensional. La Figura 14 muestra en un diagrama de bloques la estructura de control, donde 𝑓𝑐 es la conversión de la señal de control calculada en señal eléctrica (de % a 𝑉 o 𝑚𝐴) y 𝑓𝑝𝑐𝑡 es la sensibilidad del bloque conversor de la salida del transmisor (en 𝑚𝐴 o 𝑉) a un valor porcentual. Para simplificar la ilustración, se ha considerado en la figura un modelo de reactor lineal. Figura 14. Diagrama de bloques de la implementación en la GUI implementado en la GUI Nótese que 𝑓𝑐·𝐺𝑎𝑣(𝑠)=1. La referencia es introducida por el usuario en valores de ℃ que, posteriormente, es transformada en un valor porcentual por el propio controlador para obtener la señal de error. Esta estructura implica un procesamiento de los datos del sensor-transmisor y una conversión posterior de la señal de control a magnitud eléctrica, proporcionando una mayor flexibilidad al sistema de control. De esta manera, en el proceso real, sería posible emplear cualquier tipo de transmisor, incluyendo sensores con adaptadores para salida digital mediante protocolos de comunicación industriales, así como diferentes tipos de actuadores. El modelo aproximado del proceso que proporciona la herramienta de identificación paramétrica incluye la dinámica de la válvula de control, 𝐺(𝑠)=𝐺𝑣(𝑠)·𝐺𝑝(𝑠). Esta función de transferencia está aislada de los componentes eléctricos que componen el sistema, por lo que, para sintonizar un controlador PID mediante los métodos descritos en el apartado 5.1.4, debe incluirse también el producto 𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡. Teniendo en cuenta que el tiempo de respuesta del transmisor (en este caso, considerado nulo) se incluye en la función de transferencia identificada, la incorporación de este nuevo producto al modelo identificado es equivalente a realizar una transformación de (℃/%) a (%/%) en función del rango de medición del sensor mediante una sencilla multiplicación cruzada. 5.1.4 Diseño de un controlador P/PI/PID mediante reglas de sintonización Se ha diseñado una funcionalidad en la GUI para proporcionar al usuario, de forma automática, los valores del controlador P, PI o PID, en función de la opción de agresividad seleccionada, pudiendo elegir entre una respuesta agresiva, moderada o conservadora. Para ello, se ha hecho uso de reglas de sintonización basadas en ITAE [11] (para la acción proporcional) y en IMC [9] (para el control PI y PID), dadas por las siguientes expresiones:
41 Control de un reactor continuo de tanque agitado Control P 𝐾𝑐=1 2·𝐾𝑝(𝜏𝑝 𝜏𝑚𝑝)1,08 (105) Controlador PI 𝐾𝑐=𝜏𝑝 𝐾𝑝·(𝜏𝑙𝑐+𝜏𝑚𝑝) (106) 𝑇𝑖=𝜏𝑝 (107) Controlador PID 𝐾𝑐=1 𝐾𝑝(𝜏𝑝+0,5·𝜏𝑚𝑝) (𝜏𝑙𝑐+0,5·𝜏𝑚𝑝) (108) 𝑇𝑖=𝜏𝑝+0,5·𝜏𝑚𝑝 (109) 𝜏𝑝·𝜏𝑚𝑝 2𝜏𝑝+𝜏𝑚𝑝 (110) donde 𝐾𝑝 y 𝜏𝑝 son, respectivamente, la ganancia estática y constante de tiempo del sistema en lazo abierto , 𝜏𝑙𝑐 la constante de tiempo deseada en lazo cerrado y 𝜏𝑚𝑝 el retardo del proceso. Una regla comúnmente utilizada consiste en seleccionar un valor para el parámetro de ajuste de 𝜏𝑙𝑐=3𝜏 para obtener un controlador robusto y 𝜏𝑙𝑐=𝜏𝑝 para una sintonización agresiva [9]. Aunque en la literatura también existen otras reglas, en este caso, se ha optado por determinar de forma experimental tres posibilidades para 𝜏𝑙𝑐: Agresivo: 𝜏𝑙𝑐=0,5·𝜏𝑝 Moderado: 𝜏𝑙𝑐=𝜏𝑝 Conservador: 𝜏𝑙𝑐=2·𝜏𝑝 Estos valores se han seleccionado tras analizar el comportamiento de los reactores con serpentín y con camisa de refrigeración para distintos valores de la constante de tiempo en lazo cerrado deseada. El modelo lineal utilizado para calcular los parámetros de control puede modificarse desde la ventana dedicada a la identificación. La regla utilizada para el cálculo de 𝜏𝑙𝑐 no es modificable, pero se da la posibilidad al usuario de introducir manualmente los parámetros de control, pudiendo utilizarse de manera externa cualquier regla de sintonización Para proporcionar al usuario una solución, se ha diseñado un controlador PI a partir de un modelo lineal de primer orden aproximado, obtenido mediante el método de identificación descrito en el capítulo 4, tras haber realizado un ensayo en escalón abarcando todo el rango de operación del reactor. Los parámetros diseñados son mostrados al usuario por defecto en la ventana de configuración del controlador. En este caso, se han realizado varias simulaciones para distintos valores de 𝜏𝑙𝑐 próximos a 𝜏𝑝, con el objetivo de obtener un controlador a partir de un criterio de compromiso entre velocidad y suavidad en la respuesta del proceso. El retraso introducido por el mantenedor de orden cero no se ha tenido en cuenta, ya que se considera despreciable en comparación con el tiempo de respuesta del proceso.
Sistema de control por realimentación 42 Con la solución diseñada, tanto en el caso del reactor con serpentín como en el reactor encamisado, se obtienen unos márgenes de ganancia y fase de 𝑀𝐺>3,16 (10 𝑑𝐵) y 𝑀𝑓>60° para el sistema controlado en lazo abierto, dotando al control de una gran robustez. En la Figura 15, se muestra la respuesta temporal del reactor encamisado con un controlador PI ajustado mediante una sintonización IMC, para distintos valores de 𝜏𝑙𝑐. Figura 15. Comparación de desempeño entre distintos controladores PI obtenidos mediante sintonización IMC Otro método alternativo de sintonización PI es el de SIMC [12], en el que la ganancia proporcional y el tiempo integral se obtienen mediante las ecuaciones (106) y (111), respectivamente. 𝑇𝑖=min{𝜏𝑝,4(𝜏𝑙𝑐+𝜏𝑚𝑝 )} (111) donde generalmente, 𝜏𝑙𝑐≥𝜏𝑚𝑝 y 𝜏𝑙𝑐>0,1 𝜏𝑝. Éste método se ha comparado con la regla IMC para los tres niveles de agresividad definidos en la GUI, a partir de modelos lineales aproximados obtenidos por identificación paramétrica. Los resultados calculados mediante ambos métodos resultan ser idénticos, por lo que se ha optado por incluir en la GUI únicamente la regla IMC. 5.1.5 Acción anticipativa El valor de la temperatura de entrada de caudal de refrigeración puede variar, especialmente si el fluido procede de etapas anteriores en el proceso de producción. Estas variaciones son, a nivel de control, perturbaciones que afectarán a la salida del sistema y que, dependiendo del controlador implementado, podrán o no ser rechazadas por el lazo de realimentación. Para reducir los efectos sobre la señal de salida de estas perturbaciones, se ha implementado en la interfaz gráfica la opción de añadir acción anticipativa (“feedforward”) a un control PI. El esquema general de un controlador en lazo cerrado con acción anticipativa se representa mediante la Figura 16:
43 Control de un reactor continuo de tanque agitado Figura 16. Diagrama de bloques del control en lazo cerrado con acción anticipativa La ley de control viene dada en tiempo continuo por: 𝑢(𝑠)=𝐺𝑇𝐷(𝑠)· 𝐺𝐹𝐹(𝑠)·𝑑(𝑠)+𝐺𝑐(𝑠)·𝑒𝑟(𝑠) (112) donde 𝐺𝑇𝐷(𝑠) es la función de transferencia de un instrumento de medida de la perturbación 𝑑(𝑠) y 𝐺𝐹𝐹(𝑠) la función de transferencia del controlador “feedforward”. La relación 𝐺(𝑠) incluye la dinámica del proceso, la de la válvula de control y la de su actuador. La salida viene dada por: 𝑦(𝑠)=𝐺(𝑠)·𝑢(𝑠)+𝐺𝑑(𝑠)·𝑑(𝑠) (113) donde 𝐺(𝑠)=𝐺𝑎𝑣(𝑠)·𝐺𝑣(𝑠)·𝐺𝑝(𝑠). Considerando que el sistema está en equilibrio (𝑒𝑟(𝑠)=0), la salida deberá ser nula para cualquier valor que tome la señal de perturbación, por lo que, sustituyendo, se tiene: 𝐺𝐹𝐹(𝑠)=− 𝐺𝑑(𝑠) 𝐺(𝑠)·𝐺𝑇𝐷(𝑠) (114) Para una implementación más sencilla, es posible considerar únicamente la ganancia estática de la función de transferencia resultante. 𝐺𝐹𝐹(𝑠)=− 𝐾𝑑 𝐾·𝐾𝑇𝐷=− 𝐾𝑑 𝐾𝑎𝑣·𝐾𝑣·𝐾𝑝·𝐾𝑇𝐷 (115) Para el diseño implementado en la interfaz gráfica se ha usado el algoritmo de identificación paramétrica descrito en el apartado 4, con el objetivo de identificar ambas dinámicas y se ha discretizado la función de transferencia 𝐺𝐹𝐹(𝑠) resultante, mediante el método de mantenedor de orden cero (“Zero-Order Hold”), teniendo en cuenta la compensación de tiempo muerto. Además, se ha calculado la ganancia estática del controlador “feedforward”, de forma que el usuario puede elegir si implementarlo de forma estática o dinámica. En la Figura 17 se muestra una comparativa entre las respuestas temporales obtenidas para los casos de control PI con y sin acción anticipativa, aplicados a los dos modelos de reactor CSTR en el punto de operación nominal. Los resultados se han obtenido aplicando un aumento de 10 ℃ en la temperatura 𝑇𝑐𝑒 en el instante 𝑡=2 𝑚𝑖𝑛.
Sistema de control por realimentación 44 Figura 17. Comparación de las respuestas temporales ante una perturbación para los diferentes modelos controlados por un PI con y sin acción anticipativa Nótese que, al tratarse de modelos diferentes, aunque los valores iniciales de entrada (𝐹𝑐0 𝑦 𝑇𝑐𝑒0) sean iguales (véase Tablas 1 y 2), cada simulación comienza en un punto de equilibrio diferente. Desde la GUI, antes de realizar la sintonización de la acción anticipativa pueden llevarse a cabo las siguientes acciones: Modificar la función de transferencia que relaciona la entrada del proceso con la temperatura del reactor desde la ventana de “Ensayos”. Modificar la función de transferencia que relaciona la entrada de perturbación con la temperatura del reactor desde ventana de “Ensayos”. El usuario, a través de una ventana dedicada, dispone de diversas opciones para sintonizar la acción anticipativa: Indicar qué función de transferencia aproxima la relación entre la temperatura del reactor y la entrada de perturbación. Indicar qué función de transferencia aproxima la relación entre la temperatura del reactor y la variable manipulada. Cargar desde la memoria de la GUI o importar desde un archivo de texto la función de transferencia que aproxima la relación entre la temperatura del reactor y la entrada de perturbación. Cargar desde la memoria de la GUI o importar desde un archivo de texto la función de transferencia que aproxima la relación entre la temperatura del reactor y la variable manipulada. Establecer la dinámica de la acción anticipativa.
45 Control de un reactor continuo de tanque agitado Deducir automáticamente la dinámica de la acción anticipativa. Evaluar el controlador para distintos puntos de operación. Seleccionar si se desea que el algoritmo a aplicar sea dinámico o estático. La Figura 18 muestra en un diagrama de bloques la estructura de control con acción anticipativa implementada en la GUI. Esta estructura implica el procesamiento de los datos de los transmisores, proporcionando una mayor flexibilidad al sistema de control. Figura 18. Implementación en la GUI del control en lazo cerrado con acción anticipativa Con el objetivo de obtener un controlador de ganancia adimensional, la acción anticipativa considerada tiene también una ganancia adimensional. Para ello, la expresión del controlador queda: 𝐺𝐹𝐹(𝑠)=− 𝐺𝑑(𝑠) 𝑓𝑐·𝐺𝑎𝑣(𝑠)·𝐺𝑣(𝑠)·𝐺𝑝(𝑠)·𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡=− 𝐺𝑑(𝑠) 𝐺𝑇(𝑠)·𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡 (116) donde 𝐺𝑇(𝑠)=𝐺𝑣(𝑠)·𝐺𝑝(𝑠), ya que se tiene 𝑓𝑐·𝐺𝑎𝑣(𝑠)=1. Asimismo, puede utilizarse únicamente la ganancia estática del controlador obtenido como alternativa a la acción anticipativa dinámica. 5.1.6 Estrategia de control en cascada Es muy común en los procesos químicos el uso de una estrategia de control en cascada, con el objetivo de rechazar perturbaciones que actúan sobre el lazo interno, evitando que éstas lleguen a afectar al controlador principal [13]. Este esquema de control se compone de un controlador primario o maestro que se comunica con un controlador secundario o esclavo, enviándole los valores de consigna necesarios para realizar el seguimiento de la referencia principal. Para añadir más posibilidades de control dentro de la interfaz gráfica, se ha implementado una opción que permite utilizar esta estructura en el reactor encamisado. En este caso, el controlador maestro (PI o DMC) recibe la información del error de seguimiento de la variable controlada (temperatura del tanque 𝑇) y calcula el valor de la temperatura de refrigeración 𝑇𝑐 necesaria para alcanzar la referencia. El controlador esclavo (P o PI) recibe esta última señal y se encarga de manipular la válvula de control, de manera que se consiga el caudal de
Sistema de control por realimentación 52 En la Figura 24 se muestra el mismo análisis para el caso del rechazo de una perturbación de tipo escalón unitario aplicado en el instante 𝑡=1 𝑚𝑖𝑛. Puede comprobarse que con un valor de 𝛾=5 se obtiene una respuesta rápida sin sobreoscilación y una evolución suave de la abertura de la válvula de control. Figura 24. Desempeño en el rechazo de perturbaciones de un controlador “state feedback” con acción integral para distintos valores de 𝛾 5.2.4 Herramienta de sintonización de controladores de realimentación de estado Para facilitar al usuario la sintonización de controladores de realimentación del estado, se ha incluido en la GUI una herramienta que permite calcular la matriz de ganancia de realimentación mediante los métodos de asignación de polos y LQR (Figura 25). Los parámetros de ajuste disponibles para el método de asignación de polos son los polos deseados del sistema en lazo cerrado sin parte imaginaria. El usuario puede introducir directamente los valores o seleccionar una de las tres opciones sugeridas, que se detallan a continuación: Respuesta agresiva: Los polos corresponden a los del sistema en lazo abierto guardado en la memoria de la aplicación, multiplicado por un factor de 1,2 para el caso proporcional y 1,1 para el caso con acción integral. Respuesta moderada: Los polos corresponden a los del sistema en lazo abierto guardado en la memoria de la aplicación, multiplicado por un factor de 1 para el caso proporcional y 0,9 para el caso con acción integral. Respuesta conservadora: Los polos corresponden a los del sistema en lazo abierto guardado en la memoria de la aplicación, multiplicado por un factor de 0,85 para el caso proporcional y 0,7 para el caso con acción integral.
53 Control de un reactor continuo de tanque agitado Figura 25. Herramienta para la sintonización de controladores de realimentación de estado Estos valores se han seleccionado tras analizar el comportamiento del proceso para distintos valores de diseño, de forma que se ofrece al usuario tres soluciones de diferentes prestaciones, donde la diferencia radica en el tiempo de establecimiento y en el transitorio de la respuesta del proceso. El proceso de diseño realizado consiste en partir de los polos del proceso en lazo abierto y, posteriormente, multiplicarlos por un factor cercano a 1 para obtener diferentes respuestas. En caso de incluir acción integral, el factor que multiplica a los polos del sistema aproximado en lazo abierto es más pequeño con el objetivo de compensar la pérdida de estabilidad relativa. Por otro lado, los parámetros de ajuste disponibles para el método LQR son: 𝜇: Factor de ponderación de la energía de control. 𝛾: Factor de ponderación de la desviación de la integral del error respecto a la posición de equilibrio. No estará disponible en caso de configurar un controlador proporcional. Las matrices 𝑅 y 𝑄 se han diseñado con el objetivo de normalizar las variables del problema de optimización. Para ello, se ha dividido cada elemento de la diagonal de 𝑄 y 𝑅 por un factor normalizador, resultando en una penalización similar para todas las variable de decisión. De esta forma, los parámetros de ajuste 𝜇 y 𝜆 proporcionan una sintonización más simple e intuitiva del controlador. Esta normalización puede verse en la representación de las matrices 𝑄 y 𝑅 en la Figura 25.
Sistema de control por realimentación 54 Al igual que en el método de asignación de polos, el usuario puede introducir los valores de los parámetros o seleccionar una de las tres opciones predefinidas: Respuesta agresiva: Corresponde a un valor de 𝜇=0,25 para el caso proporcional y unos valores de 𝜇=2 y 𝛾=10 para el caso con acción integral. Respuesta moderada: Corresponde a un valor de 𝜇=1 para el caso proporcional y unos valores de 𝜇= 2 y 𝛾=5 para el caso con acción integral. Respuesta conservadora: Corresponde a un valor de 𝜇=5 para el caso proporcional y unos valores de 𝜇=2 y 𝛾=2 para el caso con acción integral. Estos valores se han seleccionado tras analizar la respuesta temporal del proceso para distintos valores de los parámetros de ajuste, tanto para el reactor encamisado como para el reactor enfriado por serpentín. El proceso de diseño realizado consiste en la obtención de una solución inicial para el control proporcional, 𝜇=1, para, posteriormente, aumentar y disminuir el valor del parámetro, obteniendo diferentes tiempos de establecimiento. En el caso del control con acción integral, se ha obtenido experimentalmente un valor de 𝜇 tal que no penalice de manera excesiva la energía de control para distintos valores de 𝛾. Este último parámetro se ha utilizado para obtener distintas respuestas temporales. 5.3 Control predictivo DMC 5.3.1 Algoritmo de control DMC para sistemas lineales El control predictivo basado en modelo (MPC) abarca diferentes métodos de control que implementan una aproximación práctica del controlador óptimo, utilizando para ello un modelo matemático del proceso, una acción de control que minimiza una función objetivo y una estrategia de horizonte deslizante. Entre las numerosas ventajas que proporcionan los algoritmos de control MPC, son particularmente interesantes la efectividad para tratar los sistemas multivariables y el problema del tiempo muerto, la capacidad de introducir acción anticipativa para las perturbaciones medibles y la posibilidad de controlar un gran número de procesos, tanto sencillos como de gran complejidad. Además, estos métodos resultan especialmente útiles cuando se conocen las referencias que van a aplicarse en el futuro [15]. Los controladores MPC presentan la siguiente estrategia (Figura 26): Se utiliza un modelo matemático para predecir los valores futuros de la salida del proceso 𝑦, a lo largo de un horizonte de predicción 𝑝. Esta información se calcula en cada instante en función de la información conocida del pasado y la secuencia de control calculada, que se aplicará en el futuro. La acción de control se calcula optimizando una función objetivo. Para ello, se considera un horizonte de control 𝑚 tal que 𝑚≤𝑝. Se utiliza una estrategia de horizonte deslizante: una vez calculada la secuencia óptima de acciones de control, se aplica la correspondiente al instante actual, eliminando el resto de elementos. El proceso se repite en cada instante de tiempo.
55 Control de un reactor continuo de tanque agitado Figura 26. Diagrama de funcionamiento de un controlador predictivo basado en modelo En este caso, para proporcionar al usuario herramientas de simulación y diseño para el estudio de nociones básicas sobre métodos de control MPC, se ha implementado en la interfaz gráfica un algoritmo DMC con estructura para sistemas SISO, cuyos principios se detallan a continuación. Dado el modelo matemático en escalón de un sistema: 𝑦(𝑡)=∑𝑔𝑖Δ𝑢(𝑡−𝑖) ∞ 𝑖=1 (138) donde 𝑔 es un vector que contiene los coeficientes del modelo en escalón. Nótese que el modelo tiene en cuenta el cambio en la entrada de un instante al siguiente y no el valor de la entrada en cada instante. La predicción de las salidas futuras viene dada por: 𝑦(𝑡+𝑘 | 𝑡)=∑𝑔𝑖Δ𝑢(𝑡+𝑘−𝑖) 𝑘 𝑖=1 +𝑓(𝑡+𝑘) (139) El término 𝑓(𝑡+𝑘) viene dado por: 𝑓(𝑡+𝑘)=𝑦𝑚(𝑡)+∑(𝑔𝑘+𝑖−𝑔𝑖)Δ𝑢(𝑡−𝑖) ∞ 𝑖=1 +∑(𝑔𝑑 𝑘+𝑖−𝑔𝑑 𝑖)Δ𝑑(𝑡−𝑖) ∞ 𝑖=1 + +∑𝑔𝑑 𝑖Δ𝑑(𝑡+𝑘−𝑖) 𝑘 𝑖=1 (140) donde 𝑔 corresponde a los coeficientes de la respuesta en escalón debido a la variable manipulada mientras que 𝑔𝑑 es debido a la perturbación.
Sistema de control por realimentación 56 A esta última expresión se le denomina “respuesta libre” (o no forzada) y es la señal equivalente a la salida del sistema, teniendo en cuenta todas las acciones de control tomadas en instantes anteriores a k y suponiendo incrementos nulos en la señal de control en los instantes posteriores. Si se trata de un proceso asintóticamente estable, es posible aproximar la suma de infinitos términos de las ecuaciones (138) y (140) estableciendo un horizonte de tiempo, 𝑁, lo suficientemente grande como para captar la dinámica del sistema. En la práctica, es preferible no seleccionar un valor excesivamente grande, ya que repercute en la carga computacional del algoritmo. Los horizontes asociados a los términos del proceso y a los términos de la perturbación pueden ser diferentes, al tratarse de dinámicas distintas (𝑁 y 𝑁𝑑 respectivamente). El otro término de la expresión de la predicción de salida en la ecuación (139), ∑𝑔𝑖Δ𝑢(𝑡+𝑘−𝑖) 𝑘𝑖=1 , corresponde a la respuesta forzada, que equivale a la respuesta del sistema suponiendo nulas todas las acciones de control en el pasado (instantes anteriores a k) y aplicando en los instantes presente y futuros una señal de control igual a la secuencia de control que se calculará, en este caso, resolviendo un problema de optimización. Expresando la predicción de salida en forma matricial, para varios instantes futuros, se tiene: 𝑦=𝐺𝑢+𝑓 (141) La secuencia de acciones de control 𝑢 se obtiene minimizando una función objetivo que depende del error de seguimiento de referencias y del esfuerzo de control: 𝐽=∑[𝑦(𝑡+𝑗 | 𝑡)−𝑤(𝑡+𝑗)]2+∑𝜆[Δ𝑢(𝑡+𝑗−1)]2 𝑚 𝑗=1 𝑝 𝑗=1 (142) donde 𝑝 se denomina “horizonte de predicción”, 𝑚 “horizonte de control”, 𝜆 es un parámetro de diseño utilizado para ponderar los esfuerzos de control y 𝑤 es el vector de las referencias futuras. Nótese que, para un horizonte de predicción infinito, la secuencia de control calculada en cada instante corresponde a la señal de control óptima. Teniendo en cuenta la simplificación utilizada de los horizontes de tiempo finitos, un adecuado valor para 𝑝 debe contemplar todo el transitorio de la respuesta en lazo cerrado, por lo que puede tratarse del número de muestras equivalentes al tiempo de establecimiento deseado. El horizonte de predicción es aquel en el que se evalúa la función de coste definida en la ecuación (142), por lo que la variación de este parámetro tendrá repercusión en la respuesta estimada, así como en el número de elementos de 𝑤 considerados, influyendo en la anticipación al seguimiento de referencias y en el transitorio de la respuesta del sistema. Por otro lado, el horizonte de control consiste en el número de incrementos de la señal de control considerados para controlar al proceso evaluado en el horizonte de predicción. Generalmente, un valor pequeño para 𝑚 ayuda a obtener un controlador internamente estable y, para plantas con retardo, es fundamental establecer un valor tal que 𝑚<𝑝. Un valor para el horizonte de control menor que el de predicción implica que la señal de control óptima será constante al final del horizonte de control. El problema de optimización tiene solución explícita dada por: 𝑢=(𝐺𝑇𝐺+𝜆𝐼)−1𝐺𝑇(𝑤−𝑓) (143) Con la solución calculada, se utiliza una estrategia de horizonte deslizante, aplicando únicamente el primer valor de la secuencia y repitiendo el cálculo en el instante siguiente.
57 Control de un reactor continuo de tanque agitado 5.3.2 Algoritmos para la implementación en la GUI de un control DMC para sistemas no lineales Para integrar el control DMC en la interfaz gráfica, de forma que se permita al usuario introducir los parámetros de ajuste (𝑝, 𝑚 y 𝜆), así como las condiciones de entrada de la simulación (valor de consigna, instante de aplicación de la consigna y valor e instante de la perturbación), se ha implementado el siguiente algoritmo para preparar los datos necesarios para la simulación: Las entradas son: 𝑟, 𝑡𝑟, 𝑇𝑠, Δ𝑡𝑘, las ganancias 𝐾𝑝 y 𝐾𝑑, las constantes de tiempo 𝜏𝑝 y 𝜏𝑑 y los tiempos muertos 𝜏𝑚𝑝 y 𝜏𝑚𝑑 de los sistemas de primer orden aproximados (proceso y perturbación), el tiempo de muestreo del controlador 𝑇𝐷𝑀𝐶, el número de elementos de la respuesta temporal del proceso 𝑁 y de la dinámica de perturbación 𝑁𝑑, el valor del horizonte de predicción 𝑝 y horizonte de control 𝑚 y el parámetro de ponderación de los esfuerzos de control 𝜆. Paso 1: A partir de los datos de entrada, se obtiene un vector, 𝑟𝑒𝑓 , que contiene todos los valores de consigna aplicados en cada instante de tiempo de muestreo del controlador DMC, a lo largo de toda la simulación, de forma que pueda ser interpretado por el algoritmo de control. Paso 2: Se crea un índice, 𝑖𝐷𝑀𝐶, y se le asigna el valor 𝑖𝐷𝑀𝐶=0. Este índice se utilizará para contar el número de veces que se ejecuta el algoritmo de control. Paso 3: Se construye, a partir de los datos de entrada del algoritmo, la estructura de la función de transferencia 𝐺(𝑠) que aproxima la dinámica del proceso. Paso 4: Se discretiza la función de transferencia mediante el método “Zero-Order Hold”, con un tiempo de muestreo igual al tiempo de muestreo del controlador DMC. Paso 5: Se realiza una simulación del sistema lineal en tiempo discreto y se obtienen suficientes elementos de la respuesta temporal para una entrada de tipo escalón unitario (𝑁+𝑝𝑚á𝑥 elementos), de forma que se capte la dinámica del proceso. Se crea un nuevo vector 𝑔 al que se le asignan los resultados obtenidos. Paso 6: Se repiten los pasos 3, 4 y 5 para 𝐺𝑑(𝑠), que aproxima la dinámica de perturbación, y se asignan los resultados a 𝑔𝑑. Paso 7: Se obtiene la matriz dinámica 𝐺 del algoritmo de control DMC a partir de los resultados obtenidos en el paso 5. Paso 8: Se calcula la parte “offline” de la secuencia de control óptima del algoritmo de control DMC, obteniendo la matriz 𝑀. Paso 9: Se extrae la primera fila de la matriz calculada en el paso anterior, que corresponde a la parte “offline” del primer valor de la secuencia de control y se le asigna al vector 𝐾𝐷𝑀𝐶. Paso 10: Se crea un vector nulo Δ𝑈 , que se utilizará en el algoritmo de control DMC para almacenar las acciones de control aplicadas en cada instante de tiempo. Las salidas son: El índice 𝑖𝐷𝑀𝐶 y los vectores 𝑟𝑒𝑓 , 𝑔, 𝑔𝑑, 𝐾𝐷𝑀𝐶 y Δ𝑈 . La Tabla 5 (Algoritmo 2: Preparación del control DMC) muestra de forma detallada el desarrollo matemático de este algoritmo.
Sistema de control por realimentación 58 Tabla 5. Algoritmo 2: Preparación del control DMC Entrada : 𝒓∈ℝ𝑵𝟐 , 𝒕𝒓∈ℝ𝑵𝟐 , 𝑲𝒅 ∈ℝ , 𝑲𝒑 ∈ℝ , 𝝉𝒅 ∈ℝ+ , 𝝉𝒑 ∈ℝ+, 𝝉𝒎𝒑 ∈ℝ+, 𝝉𝒎𝒅 ∈ℝ+, 𝑻𝒔∈ℝ+, 𝚫𝒕𝒌∈ℝ+, 𝑻𝑫𝑴𝑪 ∈ℝ+, 𝝀 ∈ℝ+, 𝑵 ∈ℕ , 𝑵𝒅 ∈ℕ , 𝒑 ∈ℕ , 𝒎 ∈ℕ Paso 1 : 𝑟𝑒𝑓 ∈ℝ𝒏. Obtener 𝑟𝑒𝑓 a partir de 𝑟, 𝑡𝑟, 𝑇𝑠, Δ𝑡𝑘 y 𝑇𝐷𝑀𝐶 Paso 2 : 𝑖𝐷𝑀𝐶=0 Paso 3 : 𝐺(𝑠)=𝐾𝑝 𝜏𝑝·𝑠+1𝑒−𝜏𝑚𝑝·𝑠, donde 𝑠∈ℂ Paso 4 : 𝐺(𝑧)=(1−𝑧−1) 𝒵{𝐺(𝑠) 𝑠} con tiempo de muestreo: 𝑡𝑚=𝑇𝐷𝑀𝐶, donde 𝑧∈ℂ Paso 5 : Obtener, al menos, 𝑁+𝑝𝑚á𝑥 elementos de la respuesta temporal de 𝐺(𝑧) para una entrada de tipo escalón unitario. Los resultados se asignan al vector 𝑔∈ℝ𝑁+𝑝𝑚á𝑥 Paso 6 Repetir los pasos 3, 4 y 5 para 𝐺𝑑(𝑠) y asignar los resultados a 𝑔𝑑∈ℝ𝑁𝑑+𝑝𝑚á𝑥 Paso 7 : 𝐺∈ℳ𝑝×𝑚 , tal que 𝐺=[𝑔10 … 0 𝑔2𝑔1… 0 ⋮ ⋮ ⋱ ⋮ 𝑔𝑝𝑔𝑝−1 … 𝑔𝑝−𝑚+1] Paso 8 : Calcular 𝑀=(𝐺𝑇𝐺+𝜆𝐼)−1𝐺𝑇 Paso 9 : 𝐾𝐷𝑀𝐶∈ℝ𝑝, 𝐾𝐷𝑀𝐶=[𝑀1,1 𝑀1,2 … 𝑀1,𝑝] Paso 10 : Δ𝑈 ∈ℝ𝑁+1, Δ𝑈 =0 Salida : 𝒊𝑫𝑴𝑪, 𝒓𝒆𝒇 , 𝒈, 𝒈𝒅, 𝑲 𝑫𝑴𝑪, 𝚫𝑼 Un criterio ampliamente extendido para seleccionar el tiempo de muestreo del controlador consiste en tomar un valor entre el 10% y el 25% del mínimo tiempo de respuesta deseado para el sistema en lazo cerrado [1]. En este caso, se ha seleccionado un tiempo de muestreo de 0,18 𝑚𝑖𝑛. Para simplificar el proceso de diseño por parte del usuario, este valor se mantiene fijo en la GUI. Una vez realizada la preparación del algoritmo DMC, cuando el usuario lo indique, se procederá a la simulación del sistema con los datos aportados. El cálculo de la señal de control en cada ciclo (Paso 11 del Algoritmo 1), viene representado por el siguiente algoritmo: Las entradas son: La variable controlada en el instante inicial 𝑦0 y en el instante actual 𝑦𝑘, los vectores 𝑟𝑒𝑓 , 𝑔, 𝑔𝑑 y Δ𝑈 , el último valor de la señal de control 𝑢𝑘−1, los horizontes de predicción 𝑝 y control 𝑚, el número de elementos de la respuesta en escalón del sistema lineal aproximado del proceso 𝑁 y de la dinámica de perturbación 𝑁𝑑, el tiempo equivalente al ciclo actual de simulación 𝑡, el tiempo de muestreo del controlador DMC, 𝑇𝐷𝑀𝐶, y el índice, 𝑖𝐷𝑀𝐶. Paso 1: Se obtiene el resto de la división entre 𝑡 y 𝑇𝐷𝑀𝐶. Si es nulo, se continúa en el siguiente paso, en caso contrario, se finaliza el algoritmo. Paso 2: Se aumenta en una unidad el valor del índice 𝑖𝐷𝑀𝐶. Paso 3: Se obtiene la salida medida incremental, restando a la salida actual 𝑦𝑘 el valor de la salida en el instante inicial 𝑦0.
59 Control de un reactor continuo de tanque agitado Paso 4: Se construye el vector de referencias futuras 𝑤, a partir del vector 𝑟𝑒𝑓 , del horizonte de predicción 𝑝 y del índice 𝑖𝐷𝑀𝐶. Paso 5: Se calcula la respuesta libre mediante la ecuación (140). Paso 6: Se desplazan los elementos del vector Δ𝑈 una posición a la izquierda. Paso 7: Se asigna al último elemento del vector Δ𝑈 la siguiente acción de control, calculada mediante la ecuación (143). Paso 8: Se actualiza la señal de control sumándole el último elemento del vector Δ𝑈 . Las salidas son: La actualización del vector Δ𝑈 , de la señal de control 𝑢𝑘 y del índice 𝑖𝐷𝑀𝐶. La Tabla 6 (Algoritmo 3: Cálculo de la acción de control óptima mediante DMC) muestra de forma detallada el desarrollo matemático de este algoritmo. Tabla 6. Algoritmo 3: Cálculo de la acción de control óptima mediante DMC Entrada : 𝒓𝒆𝒇 ∈ℝ𝒏, 𝑲 𝑫𝑴𝑪∈ℝ𝒑, 𝒈∈ℝ𝑵+𝒑𝒎á𝒙, 𝒈𝒅∈ℝ𝑵𝒅+𝒑𝒎á𝒙 , 𝚫𝑼 ∈ℝ𝑵+𝟏, 𝒖𝒌−𝟏∈ℝ, 𝒚𝟎∈ ℝ, 𝒚𝒌∈ℝ, 𝒕∈ℝ+, 𝑻𝑫𝑴𝑪∈ℝ+, 𝑵 ∈ℕ, 𝑵𝒅 ∈ℕ , 𝒑 ∈ℕ, 𝒎 ∈ℕ 𝒊𝑫𝑴𝑪∈ℕ Paso 1 : Si (𝑡 𝑚𝑜𝑑 𝑇𝐷𝑀𝐶)=0 ir a Paso 2. En otro caso Fin Paso 2 : 𝑖𝐷𝑀𝐶=𝑖𝐷𝑀𝐶+1 Paso 3 : 𝑦𝑚=𝑦𝑘−𝑦0 Paso 4 : 𝑤∈ℝ𝑝,𝑤=[𝑟𝑒𝑓𝑖𝐷𝑀𝐶+1 … 𝑟𝑒𝑓𝑖𝐷𝑀𝐶+𝑝]𝑇 Paso 5 : 𝑓=[𝑓(𝑡+1)… 𝑓(𝑡+𝑝)]𝑇 Paso 6 : Δ𝑈 𝑖=Δ𝑈 𝑖+1,∀ 𝑖<𝑁+1,𝑖∈ℕ+ Paso 7 : Δ𝑈 𝑁+1=𝐾𝐷𝑀𝐶·(𝑤−𝑓) Paso 8 : 𝑢𝑘=𝑢𝑘−1+Δ𝑈 𝑁+1 Salida : 𝚫𝑼 , 𝒖𝒌 , 𝒊𝑫𝑴𝑪 Nótese que, para el correcto funcionamiento del algoritmo, 𝑇𝐷𝑀𝐶 debe ser múltiplo de Δ𝑡𝑘, con el objetivo de obtener periódicamente un resto nulo en la división del paso 1. En el caso de que el usuario realice una simulación con la opción “Simulación Dinámica” seleccionada, los elementos del vector de referencias futuras del algoritmo DMC serán todos del mismo valor que la consigna indicada por el usuario. En la Figura 27 se muestran las respuestas temporales obtenidas tras realizar varias simulaciones del reactor encamisado controlado por un algoritmo DMC, con horizontes de predicción y control 𝑝=25 y 𝑚=10 y distintos valores del factor de ponderación de los esfuerzos de control 𝜆. Para valores mayores de 𝜆, se observa una menor agresividad en la respuesta transitoria a cambio de una reducción en el esfuerzo de control. Por otro lado, puede comprobarse la anticipación del controlador predictivo al seguimiento de referencias, obtenida mediante la resolución del problema de optimización utilizando un vector de referencias futuras. Nótese que el tiempo de muestreo considerado tiene efecto sobre la evolución del caudal de refrigeración, que está gobernado por la válvula de control, de tal modo que se aprecia un cierto rizado en las curvas.
Sistema de control por realimentación 60 Figura 27. Desempeño de un controlador DMC para distintos valores de 𝜆 La Figura 28 muestra una comparación entre las respuestas temporales obtenidas mediante distintas variaciones del algoritmo DMC aplicado al reactor encamisado, para una perturbación de 10℃ sobre el valor nominal de la temperatura de refrigeración a la entrada en el instante 𝑡=2 𝑚𝑖𝑛. Puede comprobarse que la adición de un modelo de perturbación mejora el comportamiento, y se obtiene un rechazo más eficaz y eficiente de la misma al utilizar el controlador en el lazo principal en una estructura de control en cascada. Figura 28. Comparación del rechazo de una perturbación para distintas variaciones del algoritmo DMC
61 Control de un reactor continuo de tanque agitado En la Figura 29 se muestran los resultados de un experimento realizado sobre el reactor encamisado, en el que se introduce, en primer lugar, una perturbación de 10 ℃ respecto al valor nominal de 𝑇𝑐𝑒 en el instante 𝑡=2 𝑚𝑖𝑛 y, posteriormente se aplica una referencia de 95 ℃ en 𝑡=10 𝑚𝑖𝑛. Los controladores se han sintonizado para proporcionar una respuesta agresiva y obtener así un buen rechazo de la perturbación y un rápido seguimiento de la referencia. Puede observarse cómo en el control DMC la señal de control está optimizada para obtener un seguimiento de las referencias futuras, a la vez que se tienen en cuenta los esfuerzos de control, de acuerdo a la función objetivo considerada. Figura 29. Comparación de la respuesta temporal obtenida mediante control PID y DMC
Control adaptativo y autoajustable 68 Figura 33. Diagrama de bloques de un controlador autoajustable Con relación a la etapa de estimación, muchos algoritmos resultan inadecuados para su implementación en línea, como es el caso del método de mínimos cuadrados descrito en el apartado 4.2. Esto es debido a que dicho estimador necesita una base de datos del proceso obtenida de manera “off-line” para calcular un único vector paramétrico. Por ello, es preferible utilizar el método de mínimos cuadrados recursivos (RLS), mediante el cual se actualizan los parámetros en cada periodo de muestreo con el objetivo de minimizar el error de estimación. Este método puede describirse mediante las siguientes ecuaciones [9]: 𝜃(𝑡)=𝜃(𝑡−1)+Γ(𝑡)𝜑(𝑡)𝜖(𝑡) (144) 𝜖(𝑡)=𝑦(𝑡)−𝜑(𝑡)𝑇𝜃(𝑡−1) (145) Γ(𝑡)=1 𝜆𝑜(Γ(𝑡−1)−Γ(𝑡−1)𝜑(𝑡)𝜑𝑇(𝑡)Γ(𝑡−1) 𝜆𝑜+𝜑(𝑡)𝑇Γ(𝑡−1)𝜑(𝑡)) (146) donde 𝜃 es el vector paramétrico que, al igual que en el algoritmo de mínimos cuadrados, contiene los parámetros a identificar, 𝜑 es el vector de regresión, que contiene los datos de entrada al algoritmo, Γ es la matriz de covarianza del error de estimación, y 𝜆𝑜 el factor de olvido, que modifica la actualización de la matriz de covarianza con el objetivo de dar un mayor peso a los datos recientes respecto a los antiguos. Este método de estimación paramétrica no requiere de una gran base de datos, porque sólo es necesario conocer los valores en el instante actual e instante previo de las variables de estimación. Además, a diferencia del método de mínimos cuadrados, no se necesita invertir matrices, lo que reduce el esfuerzo computacional del algoritmo. En caso de no disponer de información a priori de los parámetros, puede considerarse Γ0=∞𝐼 (idealmente). En el extremo opuesto, si el conocimiento a priori del vector paramétrico es perfecto, debe tomarse Γ0=0 (idealmente). Respecto al factor de olvido, éste puede variar entre 0 y 1 y su elección depende del comportamiento deseado: un valor cercano a 1 corresponde a una mayor robustez mientras que un valor cercano a 0 corresponde a una mayor velocidad de adaptación [9].
69 Control de un reactor continuo de tanque agitado La estructura del regresor y del vector paramétrico dependen del tipo de modelo que se desea identificar. En este caso, al igual que en el algoritmo implementado en la herramienta de identificación y descrito en el capítulo 4, se aproxima a un sistema de primer orden, con lo que se tiene: 𝜑𝑘=[𝑢𝑘−1−𝑛𝑑𝑒𝑙 −𝑦𝑘−1 ] (147) 𝜃=[𝑎𝑏] (148) donde pueden deducirse los parámetros del sistema equivalente en tiempo continuo a partir del método de transformación descrito en el apartado 4.3. Este procedimiento tiene la desventaja de que, al mantener un regresor nulo durante cierto tiempo, la matriz de covarianza crece exponencialmente. Para solucionar este problema, el algoritmo realiza exclusivamente la estimación paramétrica cuando hay suficiente excitación en la señal de control [9]. Para ello dispone de un sistema de detección de la excitación, que consiste en un filtrado de tipo paso alto de la referencia: 𝑦ℎ𝑝(𝑘)=𝑐ℎ𝑝𝑦ℎ𝑝(𝑘−1)+𝑘ℎ𝑝𝑐ℎ𝑝(𝑟(𝑘)−𝑟(𝑘−1)) (149) donde 𝑦ℎ𝑝 es la salida del filtro, 𝑦 la señal de salida del proceso, 𝑐ℎ𝑝 es la constante de ajuste del filtro y 𝑘ℎ𝑝 es una ganancia para escalar la señal de salida. El algoritmo de estimación puede diseñarse para ponerse en marcha cuando la señal de salida del filtro paso alto sea distinta de cero y detenerse en caso contrario. Sin embargo, es posible que pueda producirse un cambio en la referencia lo suficientemente pequeño como para no aportar excitación suficiente. Para evitar este problema, se han establecido unos determinados valores no nulos de salida del filtro para los cuales se activará o desactivará el algoritmo de estimación paramétrica. Estos valores se han determinado de manera experimental mediante la evaluación por simulación de la detección de excitación. Por otro lado, la variación de los parámetros de control se llevará a cabo una vez que la salida del filtro paso alto tenga un valor cercano a cero y el proceso de estimación recursiva se haya detenido. Esta medida se realiza para evitar que esta variación se produzca durante la respuesta transitoria del proceso. Además, con esto también se evita la aplicación del ajuste automático del controlador a partir de los primeros datos de identificación, los cuales tienen asociados una mayor incertidumbre en el error de estimación. En la Figura 34 se muestran los datos de una simulación del proceso controlado mediante un PI sin autoajuste, donde se aplica el algoritmo de identificación por mínimos cuadrados recursivos. La información a priori considerada es un vector paramétrico nulo (parámetros desconocidos) y una matriz de covarianza arbitrariamente grande. Los resultados se comparan con los parámetros obtenidos a partir de una identificación paramétrica fuera de línea, mediante el método de mínimos cuadrados, alrededor del mismo punto de operación. En la gráfica de la respuesta estimada, puede comprobarse la similitud entre la estimación (producto del regresor por el vector paramétrico calculado) y la respuesta del proceso simulado. Asimismo, también se representa la salida del filtro paso alto, que realiza la función de detección de la excitación, accionando e inhabilitando el algoritmo de identificación recursivo.
Control adaptativo y autoajustable 70 Figura 34. Resultados de una identificación paramétrica mediante mínimos cuadrados recursivos En el Anexo B se detallan los datos de simulación, así como los resultados numéricos obtenidos en el experimento de identificación paramétrica recursiva. Una vez realizado el proceso de identificación en un instante de tiempo determinado, deben utilizarse los parámetros obtenidos para realizar una sintonización automática del controlador. En este caso se ha considerado una estructura de control PI estándar y se ha empleado el método de sintonización IMC descrito en el apartado 5.1.4, en función de un valor de 𝜏𝑙𝑐 proporcionado por el usuario. Los parámetros del controlador calculados no son actualizados de forma directa, ya que se hacen pasar previamente por un filtro paso bajo para evitar efectos indeseados en la respuesta transitoria del proceso. El usuario deberá aportar también los valores iniciales de la constante proporcional y tiempo integral con los que se pondrá en marcha la simulación. A continuación, se detalla el algoritmo implementado en la interfaz gráfica para obtener en cada ciclo de simulación la actualización de los parámetros del controlador: Preparación fuera de línea: o Definir una matriz de covarianza inicial Γ, con valores arbitrariamente altos y un vector paramétrico 𝜃 inicial, preferiblemente obtenido a partir de una identificación con datos experimentales sobre el punto de operación nominal.
71 Control de un reactor continuo de tanque agitado o Diseñar los parámetros del controlador PI y asignar a los parámetros deseados los valores: 𝐾𝑐,𝑑𝑒𝑠=𝐾𝑐 y 𝑇𝑖,𝑑𝑒𝑠=𝑇𝑖. o Diseñar las constantes de los filtros paso bajo de la señal de entrada 𝑏𝑓, de la señal de salida 𝑎𝑓, de los parámetros del controlador 𝑝𝑙𝑝, y las constantes del filtro paso alto 𝑐ℎ𝑝 y 𝑘ℎ𝑝. Las entradas son: Los valores en el ciclo anterior de las variables Γ, 𝜃, 𝑢, 𝑢𝑙𝑝, 𝑦, 𝑦ℎ𝑝, 𝑦𝑙𝑝. Los valores en el ciclo actual de las variables 𝑟, 𝑢, 𝑦. Los parámetros 𝑦0, 𝐾𝑐,𝑑𝑒𝑠, 𝑇𝑖,𝑑𝑒𝑠, 𝜏𝑙𝑐, 𝑙1, 𝑙2, 𝜆𝑜, el valor del indicador 𝑅𝐿𝑆𝑓 y las constantes de los filtros 𝑎𝑓, 𝑏𝑓, 𝑐ℎ𝑝, 𝑘ℎ𝑝, 𝑝𝑙𝑝. Paso 1: Se actualiza el filtro paso bajo de la señal de control y de la señal de salida del sistema. Paso 2: Se actualiza el regresor de acuerdo a la ecuación (147). Paso 3: Se actualiza el filtro paso alto de la referencia, de acuerdo a la ecuación (149). Paso 4: Si 𝑅𝐿𝑆𝑓=0 y |𝑦ℎ𝑝𝑘|≥𝑙1, entonces 𝑅𝐿𝑆𝑓=1. Paso 5: Si 𝑅𝐿𝑆𝑓=1 y |𝑦ℎ𝑝𝑘|≤𝑙2, entonces 𝑅𝐿𝑆𝑓=0. Paso 6: Si 𝑅𝐿𝑆𝑓=1, entonces se actualiza la matriz Γ, de acuerdo a la ecuación (146) y se continúa en el paso 7. En caso contrario, continuar en el paso 11. Paso 7: Se actualiza el error de modelado mediante la ecuación (145). Paso 8: Se actualiza el vector paramétrico mediante la ecuación (144). Paso 9: Se obtienen los parámetros del modelo equivalente en tiempo continuo. Paso 10: Se realiza una sintonización IMC, de acuerdo a las ecuaciones (106) y (107), actualizando el valor de los parámetros deseados para el controlador. Paso 11: Si 𝑅𝐿𝑆𝑓=0, se actualiza el filtro paso bajo aplicado a los parámetros del controlador. Las salidas son: Las variables actualizadas Γ, 𝜃, 𝜑, K𝑐,𝑑𝑒𝑠, 𝑇𝑖,𝑑𝑒𝑠, 𝐾𝑐, 𝑇𝑖, 𝑢𝑙𝑝, 𝑦ℎ𝑝, 𝑦𝑙𝑝 y 𝑅𝐿𝑆𝑓. El desarrollo matemático de este algoritmo se detalla en la Tabla 8 (Algoritmo 4: Actualización de los parámetros de un PI autoajustable mediante RLS y sintonización IMC). La Figura 35 muestra los resultados obtenidos tras realizar una simulación del reactor CSTR controlado mediante un PI autoajustable por RLS y sintonización por IMC. Para las condiciones iniciales, se han diseñado unos parámetros de control que proporcionan una respuesta temporal lenta en la salida del proceso y un valor para 𝜏𝑙𝑐 de 1 minuto, con el objeto de aumentar la agresividad. La respuesta del proceso se hace variar entre 80 y 90 ℃ para activar el algoritmo de identificación recursiva, de forma que pueda calcularse la actualización de los parámetros de control. En el instante 𝑡=60 𝑚𝑖𝑛, se aplica una perturbación de −5 º𝐶 respecto al valor nominal de la temperatura de entrada de refrigeración. Puede comprobarse que, a medida que avanza la simulación, los valores de 𝐾𝑐 y 𝑇𝑖 evolucionan para ajustarse a la especificación de diseño, 𝜏𝑙𝑐.
Control adaptativo y autoajustable 72 Tabla 8. Algoritmo 4: Actualización de los parámetros de un PI autoajustable mediante RLS y sintonización IMC Entrada : 𝚪𝒌−𝟏∈𝓜𝒏×𝒎(ℝ), 𝜽𝒌−𝟏∈ℝ𝟐, 𝒓𝒌∈ℝ, 𝒓𝒌−𝟏∈ℝ, 𝒖𝒌∈ℝ, 𝒖𝒌−𝟏∈ℝ, 𝒖𝒍𝒑𝒌−𝟏∈ℝ, 𝒚𝒐∈ℝ, 𝒚𝒌∈ℝ, 𝒚𝒌−𝟏∈ℝ, 𝒚𝒉𝒑𝒌−𝟏∈ℝ, 𝒚𝒍𝒑𝒌−𝟏∈ℝ, 𝑲𝒄,𝒅𝒆𝒔∈ℝ, 𝝉𝒍𝒄∈ℝ, 𝒂𝒇∈ℝ+, 𝒃𝒇∈ℝ+, 𝒄𝒉𝒑∈ℝ+, 𝒌𝒉𝒑∈ℝ+, 𝒑𝒍𝒑∈ℝ+, 𝒍𝟏∈ℝ+, 𝒍𝟐∈ℝ+, 𝑻𝒊,𝒅𝒆𝒔∈ℝ+, 𝝀𝒐∈ ℝ+:𝝀𝒐∈ℤ+:𝝀𝒐∈[𝟎,𝟏], 𝑹𝑳𝑺𝒇∈ℤ+:𝑹𝑳𝑺𝒇∈[𝟎,𝟏] Paso 1 : Filtrado de la señal de entrada: 𝑢𝑙𝑝𝑘=𝑏𝑓𝑢𝑙𝑝𝑘−1+(1−𝑏𝑓)·(𝑢𝑘−𝑢0). Filtrado de la señal de salida: 𝑦𝑙𝑝𝑘=𝑎𝑓𝑦𝑙𝑝𝑘−1+(1−𝑎𝑓)·(𝑦𝑘−𝑦0). Paso 2 : Actualización del regresor: 𝜑1=−𝑦𝑙𝑝𝑘−1, 𝜑2=𝑢𝑙𝑝𝑘−1−𝑛𝑑𝑒𝑙 Paso 3 : Filtro paso alto: 𝑦ℎ𝑝𝑘=𝑐ℎ𝑝𝑦ℎ𝑝𝑘−1+𝑐ℎ𝑝𝑘ℎ𝑝(𝑟𝑘−𝑟𝑘−1) Paso 4 : Si 𝑅𝐿𝑆𝑓=0 y |𝑦ℎ𝑝𝑘|≥𝑙1, entonces 𝑅𝐿𝑆𝑓=1. Paso 5 : Si 𝑅𝐿𝑆𝑓=1 y |𝑦ℎ𝑝𝑘|≤𝑙2, entonces 𝑅𝐿𝑆𝑓=0 Paso 6 : Si 𝑅𝐿𝑆𝑓=1, entonces Γ𝑘=(1 𝜆𝑜)(Γ𝑘−1−Γ𝑘−1𝜑𝜑𝑇Γ𝑘−1 𝜆𝑜+𝜑𝑇Γ𝑘−1𝜑) y continuar en paso 7. En caso contrario, continuar en paso 11 Paso 7 : 𝜖𝑘=𝑦𝑙𝑝𝑘−𝜑𝑇𝜃𝑘−1 Paso 8 : 𝜃𝑘=𝜃𝑘−1+Γ𝑘𝜑𝜖𝑘 Paso 9 : 𝐾𝑐=𝜃2 1+𝜃1, 𝜏= Δ𝑡𝑘 |−ln(−𝜃1)|, 𝜏𝑚𝑝=𝑛𝑑𝑒𝑙·Δ𝑡𝑘 Paso 10 : 𝐾𝑐,𝑑𝑒𝑠=𝜏𝑝 𝐾𝑝(𝜏𝑙𝑐+τmp) , 𝑇𝑖,𝑑𝑒𝑠=𝜏𝑝 Paso 11 : Si 𝑅𝐿𝑆𝑓=0, 𝐾𝑐=𝑝𝑙𝑝𝐾𝑐+(1−𝑝𝑙𝑝)𝐾𝑐,𝑑𝑒𝑠, 𝑇𝑖=𝑝𝑙𝑝𝑇𝑖+(1−𝑝𝑙𝑝)𝑇𝑖,𝑑𝑒𝑠 Salida : 𝚪𝒌, 𝜽𝒌, 𝝋 , 𝒌𝒄,𝒅𝒆𝒔, 𝑻𝒊,𝒅𝒆𝒔, 𝒌𝒄, 𝑻𝒊, 𝒖𝒍𝒑𝒌, 𝒚𝒉𝒑𝒌, 𝒚𝒍𝒑𝒌, 𝑹𝑳𝑺𝒇 En la aplicación, se ha configurado un vector paramétrico inicial, el cual se ha obtenido a partir de una identificación fuera de línea con el algoritmo de mínimos cuadrados, realizando un experimento alrededor del punto de operación nominal. La matriz de covarianza inicial, así como las constantes de los filtros también han sido diseñadas experimentalmente y no son modificables por el usuario. Por otro lado, para proporcionar una solución por defecto, se han diseñado todos los parámetros de ajuste. En primer lugar, se ha sintonizado un controlador PI mediante las reglas basadas en IMC descritas en el apartado 5.1.4, partiendo del modelo obtenido con una identificación paramétrica por mínimos cuadrados. La constante del sistema en lazo cerrado puede variar en función de la agresividad deseada. En este caso, se ha considerado un valor de 1 𝑚𝑖𝑛, que da lugar a una respuesta moderada-agresiva. Respecto al factor de olvido, se ha priorizado la robustez frente a la velocidad de estimación, por lo que se ha considerado un valor de 0,995. En la aplicación, el usuario puede configurar los parámetros del controlador PI inicial, el tiempo deseado para el sistema en lazo cerrado y el factor de olvido.
73 Control de un reactor continuo de tanque agitado Figura 35. Respuesta temporal y evolución de los parámetros 𝑘𝑐 y 𝑇𝑖 para un control PI autoajustable, aplicado al caso del reactor CSTR encamisado
75 7 MEJORAS FUTURAS n este capítulo se discuten los aspectos que han quedado fuera del alcance de este trabajo, y que contribuirían a mejorar la aplicación diseñada, proporcionando nuevas herramientas de análisis y diferentes posibilidades en la configuración de la simulación y del lazo de control. Estas mejoras podrán ser incluidas en el desarrollo de una nueva versión de la interfaz gráfica, manteniendo, además, las funcionalidades ya existentes y expuestas en los capítulos anteriores. 7.1 Opciones de simulación En lo referente a las opciones de simulación, la adición de modelos de reactor de mayor complejidad permitiría la inclusión de herramientas dedicadas a un campo más avanzado de la teoría del control automático: el control multivariable. Diversos algoritmos de control podrían ser programados para facilitar al usuario el diseño de los parámetros y la observación y comparación de las respuestas obtenidas en cada caso. Otra posible mejora consiste en el escalado de las dimensiones de los reactores a valores industriales, de forma que se obtengan casos más realistas donde, además, el tiempo muerto pueda ser mayor que los valores obtenidos en este TFM. Para ello, deben modificarse los parámetros de los modelos matemáticos implementados en la aplicación, así como las unidades y límites en los campos editables de las ventanas gráficas afectados por dicha modificación. Además, debe ajustarse de nuevo cada controlador para proporcionar un valor por defecto al usuario. Para facilitar al usuario la comparación entre las diferentes respuestas temporales obtenidas con los distintos controladores, se propone introducir una opción que permita la representación, en una misma gráfica, de varias simulaciones consecutivas, mostrando, a su vez, una leyenda donde se indique a qué simulación corresponde cada curva. Respecto a la simulación dinámica, se podría incluir un panel donde se programará con antelación un conjunto de datos de entrada, definidos en distintos instantes de tiempo. De esta forma, los cambios prefijados se aplicarían automáticamente a medida que avanza la simulación, proporcionando una mayor precisión y comodidad en el ensayo. El algoritmo de control DMC, en este caso, se ejecutaría recibiendo como argumento de entrada el vector de referencias futuras, lo que permitiría una anticipación óptima de los cambios en los valores de consigna. 7.2 Identificación Para mejorar la herramienta de identificación, se propone añadir la posibilidad de seleccionar entre las distintas estructuras de modelos descritas en el apartado 4.1, así como implementar una opción que permita seleccionar dos conjuntos de datos: uno para identificación y otro para validación. Debido a esta nueva variedad de opciones, resultaría conveniente realizar un análisis estadístico basado en los residuos, mostrando al usuario en una gráfica la autocorrelación de los residuos y permitiéndole comparar los resultados con los de otro experimento diferente. Por otro lado, el algoritmo de identificación puede mejorarse considerando el número de pasos de retraso como otro elemento del vector paramétrico, lo que conlleva a modificar el problema de optimización por uno de mayor complejidad. E
Mejoras futuras 76 7.3 Opciones de control En relación a las opciones de control, se propone la adición de variaciones en el control PID, como es el caso de la estructura en paralelo u otras variantes comerciales del algoritmo, para lo cual, debería llevarse a cabo un diseño de los parámetros con la finalidad de proporcionar al usuario unos valores por defecto. Asimismo, podrían facilitarse diferentes métodos de sintonización y permitir la modificación de las reglas para los niveles de agresividad establecidos. En caso de escalar las dimensiones del reactor y obtener un modelo con un retardo de tiempo dominante (𝜏𝑚𝑝≫𝜏𝑝), resultaría conveniente añadir una estructura de predictor de Smith, con parámetros ajustables por el usuario. Por otro lado, también resultaría conveniente la posibilidad de extender la herramienta de control adaptativo y autoajustable al caso del reactor con serpentín. 7.4 Visualización Para mejorar la flexibilidad de las ventanas gráficas, se propone habilitar la opción de modificar el tamaño de las mismas, de forma que puedan expandirse o contraerse para obtener cualquier relación de aspecto, así como poder adaptarse a distintas resoluciones de pantalla. Para ello, debería diseñarse un algoritmo que calcule el tamaño inicial de los elementos y figuras presentes en cada ventana, el tamaño de fuente inicial y establecer sus posiciones relativas iniciales. Asimismo, al ejecutar la aplicación, debería determinarse la resolución del monitor y del valor del zoom seleccionado en la configuración de pantalla, de forma que se establezca una posición, tamaño y relación de aspecto inicial para todas las ventanas que componen la GUI. Por otro lado, sería necesario diseñar otro algoritmo que modifique el tamaño de los elementos, posiciones relativas y tamaño de fuente en función de la medida en píxeles del largo y ancho de cada ventana en cada instante de tiempo. Para no comprometer la claridad visual con los cambios realizados en las dimensiones, sería conveniente diseñar distintas distribuciones para los elementos de cada ventana, de manera que se adapten a distintas relaciones de aspecto, siendo de especial interés las siguientes: Relación 1:1 Relación 3:2 Relación 4:3 Relación 16:9 Actualmente, la GUI es compatible con las siguientes resoluciones: 1920x1080 (con 100% y 125% de zoom) 1400x1050 (con 100% y 125% de zoom) 1680x1050 (con 100% y 125% de zoom) 1280x1024 1280x960 1600x900 1440x900
77 REFERENCIAS [1] Mathworks, «es.mathworks.com,» [En línea]. Available: https://es.mathworks.com. [Último acceso: 2022]. [2] E. Peña, A. Pérez, A. J. Miranda y J. Sánchez, «Modelado de un reactor químico tipo CSTR y evaluación del control predictivo aplicando Matlab-Simulink,» Revista INGENIERÍA UC, vol. 15, nº 3, pp. 97-112, 2008. [3] G. L. Foutch y A. H. Johannes, «Reactors in Process Engineering,» de Encyclopedia of Physical Science and Technology, Academic Press, 2003, pp. 23-43. [4] K.-Y. Show y D.-J. Lee, «Bioreactor and Bioprocess Design,» de Biohydrogen, Elsevier, 2013, pp. 317- 337. [5] D. E. Seborg, T. F. Edgar, D. A. Mellichamp y F. J. Doyle III, Process Dynamics and Control, 4th Edition, John Wiley & Sons Ltd., 2016. [6] S. Bittanti, Model Identification and Data Analisys, John Wiley & Sons, 2019. [7] Y. Okuyama, Discrete Control Systems, Springer Science+Business Media, 2014. [8] K. Ogata, Ingeniería de Control Moderna, 5ª Edición, Pearson Education, 2010. [9] K. Johan Åström y T. Hägglund, Control PID Avanzado, Pearson Education, 2009. [10] K. J. Åström y B. Wittenmark, Computer-Controlled Systems: Theory and Design, Prentice Hall, 1997. [11] C. A. Smith y A. B. Corripio, Principles and Practice of Automatic Process Control, John Wiley and Sons, 1997. [12] S. Skogestad, «Tuning for Smooth PID Control with Acceptable Disturbance Rejection,» Industrial & Engineering Chemistry Research, vol. 45, pp. 7817-7822, 2006. [13] A. O. Ahmed, G. A. Gasmelseed, A. B. Karama y A. E. Musa, «Cascade Control of a Continuous Stirred Tank Reactor (CSTR),» Journal of Applied and Industrial Sciences, 2013. [14] A. J. Shaiju y I. R. Petersen, «Formulas for Discrete Time LQR, LQG, LEQG and Minimax LQG Optimal Control Problems,» IFAC-PapersOnLine, vol. 41, pp. 8773-8778, 2008. [15] E. Fernández Camacho y C. Bordons, Model Predictive Control, Springer Science+Business Media, 2007. [16] F. Rodríguez Rubio y M. J. López Sánchez, Control Adaptativo y Robusto, Editorial Universidad de Sevilla, 1996.
Anexo A. Manual de la interfaz gráfica de usuario 84 Señales eléctricas: Muestra la evolución temporal de las variables correspondientes a las señales eléctricas (salida de sensores y salida del controlador) (Figura A9). Figura A9. Evolución de las señales eléctricas obtenidas mediante el botón “Señales eléctricas” En la Tabla A1 se muestran los rangos de medición de los sensores-transmisores. Tabla A1. Rango de medición de los sensores-transmisores Variable Límite inferior Límite superior Temperatura del reactor 66,85 (º𝐶) 121,85 (º𝐶) Temperatura de refrigeración 66,85 (º𝐶) 121,85 (º𝐶) Concentración de reactivo a la salida 0 𝑚𝑜𝑙/𝑙 0,13 𝑚𝑜𝑙/𝑙 Temperatura de refrigeración a la entrada 16,85 (º𝐶) 61,85 (º𝐶) 3. Configuración de la simulación Para abrir el menú de configuración, haga click sobre el botón “Configuración de Simulación” en la barra superior. Se mostrará la siguiente ventana (Figura A10):
85 Control de un reactor continuo de tanque agitado Figura A10. Ventana gráfica de configuración de la simulación En el primer desplegable se puede cambiar el caso de estudio. Hay un total de 3 casos para simular: Reactor CSTR refrigerado mediante serpentín Reactor CSTR encamisado Reactor CSTR encamisado con control en cascada En el segundo desplegable se permite cambiar el método de integración numérica. Las opciones son: Método de Euler Método de Runge-Kutta de orden 4 Con el desplegable “Tipo de simulación” se pueden seleccionar dos opciones: Simulación estática: Se simula el sistema a lo largo del tiempo establecido en la casilla “Tiempo de simulación”. Los resultados se muestran en las gráficas tan pronto como finalice el bucle de simulación. Simulación dinámica: Cada ciclo de simulación y representación de datos se realiza con un tiempo de espaciamiento. Permite la simulación en tiempo real o acelerado en función del multiplicador indicado en la casilla de “Factor de tiempo”. Los dos últimos desplegables corresponden, respectivamente, a la magnitud y rango de la señal de salida de los sensores y de la señal de salida del controlador. Las posibilidades son: De 0 a 10 V De 4 a 20 mA En el panel inferior es posible modificar la desviación típica del ruido de medida simulado para la temperatura del reactor y de entrada al reactor, concentración de reactivo a la salida y temperatura de refrigeración a la entrada y a la salida de la camisa. Esta desviación típica hace referencia a las magnitudes de las señales de transmisión eléctricas (V o mA) y no a las magnitudes físicas (ºC y mol/L).
Anexo A. Manual de la interfaz gráfica de usuario 86 A la derecha, se dispone de tres campos editables donde es posible configurar con un valor numérico las siguientes opciones: Tiempo de simulación (min): Horizonte de tiempo en el que se simulan los datos en el algoritmo de simulación estática. El valor indicado debe estar entre 1 min y 200 min. Factor de tiempo: Relación entre el tiempo simulado y el tiempo real. Determina la velocidad de cálculo del simulador. El valor indicado debe estar entre 1 y 120. Al establecer un factor de tiempo de 1, los ciclos del bucle de simulación se ejecutan con un periodo igual al tiempo de muestreo, que es de 1,2 𝑠 (simulación en tiempo real). Un factor de tiempo de 120 implica reducir el periodo a 0,01 𝑠. Horizonte de visualización de datos (min): Horizonte de tiempo representado en las gráficas en el algoritmo de simulación dinámica. El valor indicado debe estar entre 1 min y 60 min. Tamaño del historial (min): Tamaño del historial donde se registran los datos en el algoritmo de simulación dinámica. El valor indicado debe estar entre 10 min y 120 min. El botón “Restaurar TF” devolverá a sus valores por defecto los parámetros de las funciones de transferencia guardadas en memoria. Para restablecer las opciones por defecto, haga click en el botón “Restablecer” y, posteriormente, acepte los cambios (esto no afectará a los modelos lineales guardados en memoria). 4. Simulación del proceso 4.1. Simulación I. Algoritmo de simulación estática Para realizar una simulación, en primer lugar, consigne los valores de entrada deseados en los campos editables, así como los instantes en los que desea aplicarlos y haga click en el botón de simulación, situado en la barra superior. Por defecto, estos campos contienen los valores correspondientes al punto de equilibrio desde el cual partirá la simulación. Posteriormente, se mostrará en las gráficas la evolución de los estados y se habilitarán los botones “Variables de proceso”, “Índices de error” y “Señales eléctricas”. Sobre las gráficas se dispone de dos botones, “Restaurar Vista” y “Ajustar Vista” (Figura A11), que pueden utilizarse para acercar o alejar el zoom automáticamente. Al situar el ratón sobre las gráficas, se despliega un menú en la parte superior (Figura A11) desde donde se permite exportarlas a formatos .png, .jpg, .tif y .pdf. También se permiten las opciones de: copiar como imagen, copiar como gráfico vectorial, desplazamiento, zoom y restaurar vista. Para limpiar las gráficas haga click en el botón “Reset” situado en la barra superior. Esta acción dejará preparada la interfaz para realizar una nueva simulación.
87 Control de un reactor continuo de tanque agitado Figura A11. Opciones de gráficas 4.2. Simulación II. Algoritmo de simulación dinámica En caso de haber seleccionado el modo de simulación dinámica, los cuadros editables de tiempo de la ventana principal quedarán inhabilitados. Se mostrarán dos casillas en la parte superior donde es posible modificar el factor de tiempo (FT) y el horizonte de visualización de datos (HVD) (véase apartado 3). Es posible seleccionar cualquier valor comprendido entre los siguientes límites: Rango disponible para el factor de tiempo: 𝐹𝑇=[1 ;120] (un valor de 𝐹𝑇=1 corresponde a una simulación en tiempo real). Rango disponible para el horizonte de visualización de datos: 𝐻𝑉𝐷=[1 ;60] 𝑚𝑖𝑛 Para llevar a cabo una simulación, deben introducirse los datos de entrada y, posteriormente, hacer click en el botón “Simulación”. Los campos editables se bloquearán y comenzarán a mostrarse los resultados en las gráficas. En este modo, se habilita un botón denominado “Autoajustar Vista”, que permite hacer un ajuste automático de los ejes y actualizarlo en cada ciclo de representación de datos (Figura A12). Para pausar temporalmente la simulación, haga click sobre el botón de “Pausa”. Las gráficas dejarán de mostrar resultados y se habilitarán de nuevo los campos editables. A partir de este momento, es posible modificar los datos de entrada que se tendrán en cuenta una vez se reanude el algoritmo. Al hacer click en “Simular”, continuará la simulación desde el punto anterior, con los nuevos datos actualizados que se hayan introducido durante la pausa. Para detener la simulación y resetear las gráficas haga click en el botón “Reset”. En este modo de simulación es posible realizar un guardado de datos en la memoria de la aplicación, con la finalidad de llevar a cabo una identificación paramétrica. Esta opción está disponible para las siguientes configuraciones: Cualquier caso de estudio con un controlador configurado en manual. Caso de estudio del reactor encamisado en cascada con el lazo interno controlado y el externo en modo manual. Para ello, se dispondrá de dos nuevos botones en la barra superior (Figura A12), de forma que, haciendo click en el botón “Grabar Ensayo” (representado por un círculo rojo), comenzarán a registrarse los resultados de la simulación. Durante la recopilación de datos, es posible realizar cambios en las condiciones de entrada a la simulación mediante el botón “Pausa” y “Simulación” y se puede cancelar la operación con el botón “Reset”, lo que producirá el borrado de todos los datos recopilados durante el proceso. Para finalizar exitosamente el guardado de datos se tienen las siguientes opciones:
Anexo A. Manual de la interfaz gráfica de usuario 88 Hacer click en el botón “Finalizar Grabación” (representado con un cuadrado rojo). Si se alcanza un tiempo de simulación de 120 minutos, la aplicación finalizará automáticamente la grabación. Si se introduce una séptima entrada en escalón, la aplicación finalizará automáticamente la grabación (la última entrada no se registrará). Con cualquiera de estas opciones, siempre que se hayan recopilado suficientes datos, se mostrará un mensaje confirmando que la operación se ha realizado con éxito. Los datos se almacenan temporalmente en la aplicación, por lo que se eliminarán al cerrar la ventana principal (para exportar los datos véase apartado 5.2). Una vez realizado el proceso, en la ventana de “Ensayos” se habilitará un nuevo botón que permite cargar los datos guardados por la aplicación. Figura A12. Guardado de datos en el modo de simulación dinámica 5. Realización de ensayos con entrada en escalón para identificación paramétrica Para realizar ensayos con entrada en escalón, se dispone de dos opciones: Ensayos para identificación usando la variable manipulada. Ensayos para identificación usando la variable de perturbación. Ambas opciones son accesibles desde la barra superior mediante sus botones correspondientes. La ventana gráfica de la herramienta es la siguiente (Figura A13):
89 Control de un reactor continuo de tanque agitado Figura A13. Ventana gráfica de identificación paramétrica El experimento realizado puede variar en función del tipo de caso de estudio y del lazo de control establecido. Dependiendo de la configuración de control, la variable de entrada puede ser: La abertura de la válvula de la línea de refrigeración, en caso de haber seleccionado inicialmente la herramienta de ensayos con la variable manipulada. La referencia de temperatura de refrigeración, en caso de haber seleccionado inicialmente la herramienta de ensayos con la variable manipulada y tener configurado un control en cascada con el lazo interno cerrado. La temperatura de refrigeración de entrada, en caso de haber seleccionado inicialmente la herramienta de ensayos con la variable de perturbación. La temperatura de entrada al reactor, en caso de haber seleccionado inicialmente la herramienta de ensayos con la variable de perturbación y tener configurado un control en cascada con el lazo interno cerrado. Asimismo, la variable indicadora de punto de funcionamiento puede ser: La abertura inicial de la válvula de refrigeración, en caso de control regulatorio básico y en caso de control en cascada en modo manual. La referencia inicial de temperatura de refrigeración a la salida de la camisa, en caso de control en cascada con el lazo interno cerrado.
Anexo A. Manual de la interfaz gráfica de usuario 90 5.1. Obtención de datos para identificación En el panel superior puede seleccionarse el número de entradas que se aplicarán en la simulación. A medida que se vayan agregando, aparecerán en el panel nuevos campos editables para rellenarlos con los datos del ensayo. Se dispone de dos columnas de datos de entrada, siendo la primera la correspondiente a los valores de la variable de entrada y la segunda, los instantes en los que se van a aplicar. Debajo del panel, se muestran dos campos editables en los que puede especificarse el tiempo de simulación (Ts) y un valor inicial para la variable indicadora de punto inicial de funcionamiento. Una vez configuradas las condiciones del ensayo, haga click en el botón simular. Los resultados se mostrarán en las gráficas de la derecha. Otra forma de obtener datos para la identificación consiste en realizar una grabación desde la ventana principal, activando el modo de simulación dinámica. Para ello, seleccione el valor inicial de la válvula de refrigeración, haga click en el botón “Guardar datos” localizado en la barra superior de la ventana principal e inicie una simulación (véase apartado 4.2). Los datos registrados previamente a la introducción del primer escalón no se utilizarán en el algoritmo de identificación. De esta manera, es posible seleccionar el punto de equilibrio en torno al cual se desea llevar a cabo el experimento. En caso de haber datos registrados en memoria, al abrir la ventana de “Ensayos” se mostrará el botón “Cargar Datos” y, haciendo click sobre él, se representarán en las gráficas los resultados recopilados durante la grabación. Los datos se cargan desde la memoria interna de la aplicación y pueden ser exportados a una hoja de cálculo o a un documento de texto (ver apartado 5.2). Asimismo, si se dispone de un archivo con datos guardados, también es posible cargarlos en la herramienta mediante el botón “Importar Datos” (véase Figura A13). Para ello, debe seleccionar el tipo de extensión (.xlsx o .txt) mediante el menú desplegable de la nueva ventana emergente (Figura A14) y, posteriormente, seleccionar el archivo deseado y hacer click en “Abrir”. Figura A14. Importación de datos Para obtener un conjunto de datos válido para la identificación, se recomienda: Introducir el primer escalón una vez el sistema se haya estabilizado. Evitar que la primera entrada sea igual al valor nominal (escalón de valor nulo).
91 Control de un reactor continuo de tanque agitado Dejar un margen de separación temporal entre distintos escalones, de forma que el sistema siempre se estabilice. Evitar que dos entradas consecutivas tengan un valor muy próximo. Haciendo click en el botón “Variables de Proceso” se generarán tres gráficas: En una se representa la evolución temporal de los estados, tanto del proceso como del modelo identificado. En otra se representa la evolución de la variable de perturbación, la potencia calorífica disipada y la constante cinética. En la última, se representan los datos del caudal de refrigeración y de porcentaje de abertura de la válvula. 5.2. Exportación de datos de simulación Tras obtener los datos con cualquiera de los métodos disponibles (simulación desde la herramienta, carga de datos o importación de datos), en la parte inferior se habilitará una nueva zona, donde se permite exportar los datos de la simulación. El usuario puede elegir el tipo de extensión del archivo, disponiendo de las siguientes opciones: Hoja de cálculo “.xlsx”. Archivo de texto “.txt” (con datos separados por comas). El nombre del archivo deberá indicarse en el correspondiente campo editable antes de llevar a cabo la exportación de datos. Al hacer click en el botón “Exportar” se abrirá una nueva ventana en la que deberá seleccionar el directorio de destino y, pulsando en “Seleccionar carpeta”, se mostrará un mensaje confirmando la ruta donde se han almacenado los datos (Figura A15). Figura A15. Mensaje de confirmación de exportación correcta de datos y directorio de ubicación 5.3. Identificación paramétrica En la parte inferior, se dispone de un panel donde es posible realizar la identificación paramétrica de una función de transferencia aproximada, a partir de los datos de la simulación. Existen dos posibilidades: Identificación única: Se identifica una única función de transferencia para cada estado interno del proceso, a partir de los datos de simulación. Identificación por puntos: Se identifica una función de transferencia por cada escalón introducido en el panel de ensayo, a partir de los datos de simulación y para cada estado interno del proceso.
Anexo A. Manual de la interfaz gráfica de usuario 92 Al hacer click sobre el botón “Identificación”, se realizará la identificación con la opción seleccionada, seguida de una simulación con la función de transferencia calculada. Los nuevos datos se mostrarán superpuestos en las gráficas (Figura A16), permitiendo observar las diferencias entre la respuesta temporal del modelo del reactor y del modelo calculado. En el algoritmo de identificación, los datos se registran a partir del primer escalón introducido, por lo que es recomendable aplicar dicho escalón en un instante que no comprometa la estabilización inicial. Por ello, en la herramienta de ensayos, el menor instante de tiempo admisible es de un minuto, aunque se recomienda que sea mayor en caso de que el valor inicial de la variable manipulada sea diferente al nominal. Asimismo, en caso de seleccionar el método de identificación por puntos, los escalones no pueden ser nulos y debe seleccionarse un primer valor para la variable manipulada diferente al valor inicial. Figura A16. Resultados de la identificación paramétrica mediante la opción “Identificación por puntos” Con el botón “Obtener TF” emergerá una nueva ventana (Figura A17) donde se mostrarán las expresiones matemáticas de las funciones de transferencia en tiempo continuo que se han obtenido en la identificación. Estas funciones engloban la dinámica del proceso (asociada a cada una de las tres salidas consideradas) y de la válvula localizada en la línea de refrigeración. En caso de haber seleccionado la opción “Identificación por puntos”, se mostrarán de forma numerada todos los modelos obtenidos a partir de los datos registrados para cada entrada en escalón.
93 Control de un reactor continuo de tanque agitado Figura A17. Modelos lineales obtenidos con una identificación paramétrica mediante la opción “Identificación por puntos” La numeración de los modelos lineales obtenidos mediante el método de identificación por puntos hace referencia al orden en el cual se han aplicado los cambios en los valores de la variable de entrada, siendo el número 1 el relativo al primer escalón introducido y el número 6, al último. En la parte superior, se indica qué relación de entrada/salida corresponde a cada columna representada, que dependerá del tipo de configuración establecida para el controlador, así como el tipo de modelo seleccionado. Con el botón “Reset”, se eliminarán los datos de la identificación, permitiendo realizarla de nuevo con cualquiera de las dos opciones disponibles. 5.4. Herramienta de análisis, guardado en memoria y exportación de funciones de transferencia En la parte inferior de la ventana, se dispone del botón “TF Manager”, mediante el cual se accede a una herramienta para analizar funciones de transferencia (Figura A18).
Anexo A. Manual de la interfaz gráfica de usuario 100 Figura A25. Configuración del tipo de acción derivativa con filtro Mediante el panel “Autoajuste”, es posible obtener soluciones basadas en reglas de sintonización IMC (o basado en ITAE en caso de sólo proporcional). La solución proporcionada dependerá de las acciones incluidas mediante los marcadores y del tipo de respuesta deseada, pudiendo elegir entre: Respuesta agresiva Respuesta moderada Respuesta conservadora Haciendo click en el botón “Bode” se mostrarán en dos nuevas ventanas el diagrama de Bode del sistema controlado en lazo abierto con los valores de los márgenes de fase y ganancia. En una ventana se representará la amplitud en valor absoluto (Figura A26) y en la otra, en decibelios (Figura 27). Figura A26. Diagrama de bode del proceso controlado en lazo abierto (I)
101 Control de un reactor continuo de tanque agitado Figura A27. Diagrama de bode del proceso controlado en lazo abierto (II) Para sintonizar manualmente un controlador basado en PID desde la ventana de configuración, con el fin de obtener una ganancia adimensional, debe utilizarse el siguiente modelo: 𝑓𝑐·𝐺𝑎𝑣(𝑠)·𝐺𝑣(𝑠)·𝐺𝑝(𝑠)·𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡 donde 𝐺𝑣(𝑠)·𝐺𝑝(𝑠) corresponde a la dinámica del proceso, incluyendo la válvula de control, y equivale al sistema obtenido en la identificación paramétrica. 𝐺𝑇𝐷(𝑠) es la relación del sensor-transmisor, que depende del tipo de señal configurada, y el factor 𝑓𝑝𝑐𝑡 es la sensibilidad del bloque conversor de la salida del transmisor (en 𝑚𝐴 o 𝑉) a un valor porcentual. El producto 𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡 puede obtenerse utilizando la Tabla A1. El controlador simulado está configurado para transformar la señal de control a la magnitud eléctrica necesaria para accionar la válvula de control (V o mA, según se tenga configurado), por lo que no es necesario incluir la dinámica del actuador de la válvula en el diseño del controlador 𝑓𝑐·𝐺𝑎𝑣(𝑠)=1. 6.2. Sintonización de la acción anticipativa En caso de seleccionar desde el menú desplegable el algoritmo PI con acción anticipativa, se mostrará la ventana representada en la Figura A28. Se dispone de dos campos editables para introducir los valores de la constante proporcional, 𝐾𝑐 y el tiempo integral 𝑇𝑖. Mediante el botón “Feedforward”, se accederá a una herramienta para configurar la acción anticipativa (Figura A29). Desde esta herramienta, en la zona izquierda de la ventana, es posible indicar las funciones de transferencia que van a tenerse en cuenta en el cálculo del controlador (relativas al proceso y al efecto de la perturbación). Para ello, rellene los campos editables o haga click en el botón “Cargar”, con lo que se rellenará automáticamente, a partir de la información guardada en la memoria de la aplicación (para modificar los modelos aproximados guardados en memoria, véase apartado 5.4). De manera análoga, es posible cargar los datos guardados en un archivo de texto (.txt) mediante el botón “Importar”. La función de transferencia del proceso indicada incluye la dinámica de la válvula de control.
Anexo A. Manual de la interfaz gráfica de usuario 102 Figura A28. Sintonización de un controlador PI con acción anticipativa A la derecha, se muestra la expresión del controlador. La sintonización es posible mediante las siguientes opciones: Introduciendo manualmente los parámetros del “feedforward” Haciendo click en el botón “Deducir del modelo”, con lo que se realizará el cálculo a partir de las funciones de transferencia indicadas en la zona izquierda de la ventana. Más abajo, se mostrará en un recuadro el valor de la ganancia estática correspondiente al controlador introducido (este campo no es modificable). Si desea llevar a cabo el proceso completo de identificación y posterior sintonización, deberá realizar los siguientes pasos: 1) En la ventana de ensayos con cambio en la variable manipulada, acceder a la herramienta de análisis de las funciones de transferencia. 2) Guardar en memoria interna o exportar datos de la función de transferencia aproximada. 3) Repetir los pasos 1 y 2 en la ventana de ensayos con cambio en la variable de perturbación. 4) Acceder a la configuración del controlador desde la ventana principal y seleccionar un controlador PI con acción anticipativa. 5) Hacer click en el botón “Feedforward”. 6) Cargar o importar el modelo asociado a la dinámica de perturbación. 7) Cargar o importar el modelo asociado a la dinámica del proceso incluyendo la correspondiente a la válvula de control. 8) Hacer click en el botón “Deducir del modelo”.
103 Control de un reactor continuo de tanque agitado Figura A29. Herramienta para sintonizar la acción anticipativa El cálculo de la acción anticipativa que se lleva a cabo es el siguiente: 𝐺𝐹𝐹(𝑠)=− 𝐺𝑑(𝑠) 𝐺𝑣(𝑠)·𝐺𝑝(𝑠)·𝐺𝑇𝐷(𝑠)·𝑓𝑝𝑐𝑡 Desde el panel “Tipo de acción anticipativa” puede seleccionarse si el controlador se configurará en modo estático o en dinámico. Modo estático: Se utiliza el dato de la ganancia estática de la acción anticipativa. Modo dinámico: Se utilizan los datos dinámicos indicados en la función de transferencia de la acción anticipativa. El panel “Previsualización” contiene los elementos necesarios para obtener una previsualización del modelo controlado, para una entrada unitaria de perturbación. Además, se dispone de un campo editable donde es posible modificar el valor de la referencia para evaluar la influencia de la perturbación en distintos puntos de funcionamiento y un marcador para activar o desactivar la acción anticipativa. Los parámetros del controlador PI utilizados en esta simulación de prueba, son los correspondientes a los campos editables de la ventana de configuración de control, (𝐾𝑐 y 𝑇𝑖), que deberán seleccionarse antes de acceder a esta herramienta. Una vez seleccionado el tipo de acción anticipativa deseada desde el panel correspondiente y haciendo click sobre el botón “Aplicar”, la ventana se cerrará y se guardarán los cambios de manera provisional en la memoria de la aplicación. Para hacer efectivos dichos cambios, deberá de aceptarse la configuración del PI desde la ventana de configuración de control, haciendo click en el botón “Aceptar”. En la ventana principal se mostrará el diagrama de proceso del reactor controlado con acción anticipativa (Figura A30).
Anexo A. Manual de la interfaz gráfica de usuario 104 Figura A30. Lazo de control configurado con acción anticipativa 6.3. Sintonización de un controlador por realimentación del estado En caso de seleccionar un algoritmo de control por realimentación del estado desde el menú desplegable, se mostrará la siguiente interfaz (Figura A31). Figura A31. Sintonización de un controlador por realimentación del estado con acción integral
105 Control de un reactor continuo de tanque agitado Los campos editables disponibles representan la ganancia aplicada a cada variable interna del proceso y a la integral del error (en caso de incluir acción integral). Mediante el botón “Diseño” se accede a una herramienta para diseñar los parámetros del controlador mediante los métodos de diseño de asignación de polos y LQR (Figura A32). Una vez se haya accedido a la herramienta, deberá seleccionarse un método de diseño desde el panel “Método de diseño”. En función de la opción escogida, se habilitarán los campos editables correspondientes para ajustar el controlador. Para la sintonización mediante asignación de polos, deberá indicar en los campos editables correspondientes los polos deseados para el sistema en lazo cerrado. En caso de introducir dos polos iguales, no se permitirá el cálculo de los parámetros del controlador y, por lo tanto, no se aplicará ningún cambio. Mediante el panel de autoajuste se permite seleccionar entre tres tipos de soluciones predeterminadas: agresiva, moderada y conservadora. De igual forma, para sintonizar el controlador mediante LQR, deberá introducir el valor del factor de ponderación de la energía de control. Si el algoritmo seleccionado incluye acción integral, también podrá modificarse el factor de ponderación de los incrementos en la integral del error. La agresividad del controlador se verá potenciada por: Un valor elevado del parámetro 𝛾 Un valor reducido del parámetro 𝜇 Mediante el panel de previsualización, se permite obtener una simulación del proceso sin ruido para una referencia unitaria, a partir de los parámetros indicados para el controlador. Figura A32. Herramienta para sintonización de un controlador por realimentación del estado
Anexo A. Manual de la interfaz gráfica de usuario 106 Una vez seleccionado el método de diseño y sintonizado el controlador, al hacer click sobre el botón “Aplicar” se cerrará la herramienta, se calcularán los elementos de la matriz de realimentación del estado (incluyendo 𝐾𝑖 en caso de seleccionar acción integral) y se modificarán los campos editables correspondientes a los parámetros del controlador. 6.4. Sintonización de un controlador predictivo DMC Al seleccionar desde el menú desplegable el algoritmo de control predictivo DMC, se mostrará la siguiente interfaz (Figura A33): Figura A33. Sintonización de un controlador predictivo DMC Los parámetros modificables son: El horizonte de predicción “𝑝”, expresado en múltiplos del tiempo de muestreo. El horizonte de control “𝑚”, expresado en número de elementos de la secuencia de control calculada. El factor de ponderación de los esfuerzos de control “𝜆” (adimensional). El tiempo de muestreo del controlador puede observarse haciendo click en el botón de información y está fijado en un valor de 0,18 minutos. La agresividad viene determinada, principalmente, por el valor del factor “𝜆”, de forma que un valor bajo corresponde a una respuesta más rápida y un valor alto a una respuesta con menor gasto energético. Para deducir el modelo de respuesta en escalón utilizado por el algoritmo, se parte de la función de transferencia cargada en la memoria de la aplicación. Dicha función puede modificarse desde la herramienta de ensayos con entrada en escalón, accediendo a la ventana de análisis de funciones de transferencia y guardando en memoria el modelo deseado (Véase apartado 5.4).
107 Control de un reactor continuo de tanque agitado 6.5. Configuración del control en cascada En caso de haber seleccionado, desde el menú de configuración de simulación, el modelo de reactor encamisado con estructura de control en cascada, en la ventana de opciones de control se mostrarán dos desplegables en lugar de uno (Figura A34), que corresponden a los algoritmos del controlador maestro y esclavo. En este caso las opciones disponibles para el controlador maestro son: Modo manual Control P Control PI Control PI con acción anticipativa Control predictivo DMC Para el controlador esclavo, las opciones disponibles son el modo manual, el control P y el control PI. Figura A34. Ventana gráfica de opciones de control para el reactor encamisado con estructura de control en cascada Si se configura un control P o PI para el lazo interno y el controlador primario se configura en manual, el sistema realizará un seguimiento de referencias para la variable de temperatura de refrigeración. En el caso opuesto, seleccionando un algoritmo para el primario y configurando el modo manual en el secundario, el sistema se configurará en lazo abierto, permitiendo manipular la válvula de refrigeración. Para llevar a cabo una sintonización por IMC/ITAE a partir de modelos previamente identificados, deben realizarse los siguientes pasos: 1) Configurar el controlador esclavo en modo manual. 2) Acceder a la herramienta de ensayos con cambios en la variable manipulada. 3) Realizar una identificación (opcional).
Anexo A. Manual de la interfaz gráfica de usuario 108 4) Acceder a la ventana de análisis de funciones de transferencia. 5) Guardar en la memoria interna la función de transferencia aproximada. 6) Acceder a la ventana de configuración de control desde la ventana principal. 7) Configurar un controlador esclavo. En caso de seleccionar el control PI, éste podrá sintonizarse automáticamente a partir del modelo guardado anteriormente (Botón “IMC”). Si, por el contrario, el algoritmo seleccionado es el control P, los pasos anteriores serán opcionales, ya que deberá ajustarse manualmente, para lo cual, se recomienda sintonizar una ganancia elevada, pero que estabilice al sistema para referencias de temperatura de refrigeración. 8) Acceder de nuevo a la herramienta de ensayos con cambios en la variable manipulada, que ahora será la referencia de la temperatura de refrigeración. 9) Repetir los pasos 3 (opcional), 4, 5 y 6. 10) Configurar un controlador maestro. Este controlador podrá sintonizarse automáticamente a partir del modelo guardado anteriormente. En caso de seleccionar un controlador DMC, no se dispondrá de la opción de sintonización automática y el modelo de predicción utilizado será el obtenido en el paso 9. De manera análoga, partiendo del paso 7, es posible sintonizar la acción anticipativa para la variable de perturbación externa: 8) Acceder a la herramienta de ensayos con cambios en la variable de perturbación, que ahora será la temperatura de entrada al reactor. 9) Repetir los pasos 3 (opcional), 4, 5 y 6. 10) Seleccionar la opción “Control PI con feedforward” para el controlador maestro, y hacer click en el botón “Feedforward”. 11) Cargar o importar el modelo de perturbación y el del proceso (incluyendo la dinámica de la válvula). 12) Hacer click en el botón “Deducir del modelo” 13) Hacer click en “Aplicar”. La ventana se cerrará y se guardarán los cambios temporalmente. 14) Hacer click en “Aceptar” para hacer efectivo los cambios Si no se ha llevado a cabo un guardado en memoria de las funciones de transferencia necesarias, se utilizarán las configuradas por defecto al inicio de la aplicación. El algoritmo de control DMC también tiene en cuenta un modelo de perturbación de la dinámica externa (temperatura del reactor a la entrada), que puede ser modificado desde la ventana de ensayos con cambios en la variable de perturbación (véase apartado 5.4). Haciendo click sobre el botón “Bode” se mostrará en una figura el diagrama de Bode del sistema controlado en lazo abierto. En caso de seleccionar un control P/PI para el control maestro, el diagrama de Bode se calculará cerrando el lazo interno con los parámetros indicados y dejando abierto el externo. 7. Ensayos con control adaptativo y autoajustable En caso de tener configurado el modelo encamisado, podrá accederse a una ventana destinada a realizar ensayos con controladores adaptativos y autoajustables (Figura A35), haciendo click en el botón correspondiente localizado en la barra superior.
109 Control de un reactor continuo de tanque agitado Una vez se haya accedido a la herramienta, en la parte superior, se dispone de un menú desplegable donde es posible seleccionar el algoritmo de control, existiendo tres posibilidades: Gain Scheduling PI, Gain Scheduling DMC y Controlador Autoajustable RLS. En el panel “Parámetros” se muestran los campos editables de los parámetros del controlador, que variarán en función del algoritmo de control seleccionado. En el panel inferior pueden seleccionarse los valores de entrada del ensayo (referencias de temperatura, variable de perturbación e instantes de aplicación), así como establecer el tiempo de simulación. Al hacer click en el botón “Simulación” se mostrará la evolución de los estados del proceso en las gráficas de la parte derecha de la ventana. Para observar la evolución temporal de las variables del proceso, haga click en el botón “Variables de proceso”, localizado en la parte inferior de la ventana. Mediante el botón “Variación de parámetros” se mostrará la evolución temporal de los parámetros del controlador seleccionado. Figura A35. Ventana gráfica de control adaptativo y autoajustable 7.1. Configuración del control Gain Scheduling PI Con el algoritmo Gain Scheduling PI seleccionado, se permite establecer fuera de línea un conjunto de tres controladores PI, definiendo una tabla donde se relaciona cada par de parámetros (𝐾𝑐 y 𝑇𝑖) con la zona que van a controlar. Cada una de estas zonas está determinada por un rango de temperatura, de forma que, en función de la medida de la temperatura del reactor, actuará el controlador correspondiente, de acuerdo a la información descrita en la tabla.
Anexo B. Resultados experimentales de algoritmos de identificación 116 2.2. Configuración de la simulación La configuración de la simulación, así como las condiciones iniciales para la identificación recursiva en línea se muestran en la Tabla B3. Tabla B3. Configuración de la simulación en el experimento de identificación paramétrica RLS VARIABLE DESCRIPCIÓN VALOR 𝑴𝑰 Método de integración numérica: 𝑀𝐼=1: Runge-Kutta orden 4 𝑀𝐼=2: Euler 1 𝑻𝟎 Temperatura inicial del reactor (en ºC) 87,81 𝒄𝑨𝟎 Concentración inicial de reactivo a la salida (en mol/L) 0,04 𝑻𝒄𝟎 Temperatura inicial de refrigeración (en ºC) 82 𝒒𝒄𝟎 Caudal inicial de refrigeración (L/min) 90 𝒗𝒂𝒍𝒗𝟎 Abertura inicial de la válvula de refrigeración (en %) 64,29 𝑲𝒑 Constante proporcional del controlador PI (en ΔºC/Δ%) −1,5 𝑻𝒊 Tiempo integral del controlador PI (en minutos) 0,89 𝒕𝒉𝒆𝒕𝒂𝟎 Vector paramétrico inicial [0;0] 𝑷 Matriz de covarianza inicial [100 104 104100] 𝒐𝒍𝒗𝒊𝒅𝒐 Factor de olvido (adimensional) 0,99 𝒏 Simulación con ruido: 𝑛=0: Simulación sin ruido 𝑛=1: Simulación con ruido 0 𝑴𝒐𝒅𝒆.𝑺𝑻 Tipo de sensor: 𝑀𝑜𝑑𝑒.𝑆𝑇=1: Señal de [0−10] V 𝑀𝑜𝑑𝑒.𝑆𝑇=2: Señal de [4−20] mA 1 𝑴𝒐𝒅𝒆.𝑪𝒐𝒖𝒕 Señal de salida del controlador: 𝑀𝑜𝑑𝑒.𝐶𝑜𝑢𝑡=1: Señal de [0−10] V 𝑀𝑜𝑑𝑒.𝐶𝑜𝑢𝑡=2: Señal de [4−20] mA 1 𝒕𝒓𝒆𝒇 Vector de instantes de aplicación de las referencias (en minutos) [1; 10; 20; 30; 40] 𝒓𝒆𝒇 Vector de referencias (en ºC) [92,81; 82,81; 92,81; 355,96; 92,81] 𝑻𝒔 Tiempo de simulación (en minutos) 50
117 Control de un reactor continuo de tanque agitado La configuración de la simulación y condiciones iniciales para la identificación fuera de línea se muestran en la Tabla B4. Tabla B4. Configuración de la simulación en el experimento de identificación paramétrica mediante LS (II) VARIABLE DESCRIPCIÓN VALOR 𝑴𝑰 Método de integración numérica: 𝑀𝐼=1: Runge-Kutta orden 4 𝑀𝐼=2: Euler 1 𝑻𝟎 Temperatura inicial del reactor (en ºC) 87,81 𝒄𝑨𝟎 Concentración inicial de reactivo a la salida (en mol/L) 0,042 𝑻𝒄𝟎 Temperatura inicial de refrigeración (en ºC) 82 𝒒𝒄𝟎 Caudal inicial de refrigeración (L/min) 90 𝒗𝒂𝒍𝒗𝟎 Apertura inicial de la válvula de refrigeración (en %) 64,29 𝒏 Simulación con ruido: 𝑛=0: Simulación sin ruido 𝑛=1: Simulación con ruido 0 𝑴𝒐𝒅𝒆.𝑺𝑻 Tipo de sensor: 𝑀𝑜𝑑𝑒.𝑆𝑇=1: Señal de [0−10] V 𝑀𝑜𝑑𝑒.𝑆𝑇=2: Señal de [4−20] mA 1 𝑴𝒐𝒅𝒆.𝑪𝒐𝒖𝒕 Señal de salida del controlador: 𝑀𝑜𝑑𝑒.𝐶𝑜𝑢𝑡=1: Señal de [0−10] V 𝑀𝑜𝑑𝑒.𝐶𝑜𝑢𝑡=2: Señal de [4−20] mA 1 𝒕𝒔𝒕𝒑 Vector de instantes de aplicación de cambios en la variable manipulable (en minutos) [1; 10; 20; 30; 40] 𝒔𝒕𝒑 Vector de valores de la variable manipulable (en %) [74,29; 54,29; 74,29; 54,29; 74,29] 𝑻𝒔 Tiempo de simulación (en minutos) 50
Anexo B. Resultados experimentales de algoritmos de identificación 118 2.3. Resultados Los resultados obtenidos están recogidos en la Tabla B5 y se discuten en el apartado 6.2. Tabla B5. Resultados de la identificación paramétrica RLS VARIABLE DESCRIPCIÓN VALOR 𝒂𝒎_𝑹𝑳𝑺 Valor medio del parámetro 1 obtenido mediante el método RLS, una vez eliminados los datos de estimación inicial (adimensional) −0,98 𝒃𝒎_𝑹𝑳𝑺 Valor medio del parámetro 2 obtenido mediante el método RLS, una vez eliminados los datos de estimación inicial (adimensional) −0,01 𝑲𝒎_𝑹𝑳𝑺 Valor medio de la ganancia estacionaria obtenida mediante el método RLS, una vez eliminados los datos de estimación inicial (en ΔºC/Δ%) −0,55 𝒕𝒂𝒖𝒎_𝑹𝑳𝑺 Valor medio de la constante de tiempo obtenida mediante el método RLS, una vez eliminados los datos de estimación inicial (minutos) 0,92 𝒆𝒎_𝑹𝑳𝑺 Promedio del valor absoluto del error de estimación mediante RLS, una vez eliminados los datos de estimación inicial (en ΔºC) 5·10−3 𝒂_𝑳𝑺 Valor del parámetro 1 obtenido mediante LS (adimensional) −0,98 𝒃_𝑳𝑺 Valor del parámetro 2 obtenido mediante LS (adimensional) −0,01 𝑲_𝑳𝑺 Valor de la ganancia estacionaria obtenida mediante el método LS (en ΔºC/Δ%) −0,54 𝒕𝒂𝒖_𝑳𝑺 Valor de la constante de tiempo obtenida mediante el método LS (minutos) 1,08 𝒆𝒎_𝑳𝑺 Promedio del valor absoluto del error de estimación mediante LS (en ΔºC) 0,01
119 Control de un reactor continuo de tanque agitado La Figura B2 ilustra los resultados del experimento. Figura B2. Resultados de la identificación paramétrica RLS y comparación con el método LS