Implementación de controladores predictivos en LabVIEW

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Implementación de controladores predictivos en LabVIEW Autor: José Ramón Romero Gómez Tutor: Daniel Limón Marruedo Dep. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 iii Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Implementación de controladores predictivos en LabVIEW Autor: José Ramón Romero Gómez Tutor: Daniel Limón Marruedo Profesor Titular Dep. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 v Trabajo Fin de Grado: Implementación de controladores predictivos en LabVIEW Autor: José Ramón Romero Gómez Tutor: Daniel Limón Marruedo El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2016 El Secretario del Tribunal vii A mi familia ix Agradecimientos Este proyecto representa la culminación del esfuerzo y el trabajo realizado durante los años de estudio de este Grado. En primer lugar quiero dar las gracias a mis padres, no solo por haberme dado esta la oportunidad, sino por su incondicional apoyo. De igual forma quiero agradecer a toma mi familia su apoyo y confianza en mí para la realización de estos estudios. Agradecer también a mis amigos y a mis compañeros de clase haber sido una gran ayuda durante el transcurso de estos años. Quiero agradecer personalmente a mi tutor Daniel Limón Marruedo por su tiempo, su dedicación y su cercanía. Todo ello ha hecho posible la realización de este Trabajo Fin de Grado. De igual forma, agradecer a Pablo Krupa su tiempo y ayuda en los inicios de este proyecto. José Ramón Romero Gómez Sevilla, 2016 2.3 Control Predictivo Basado en Modelo que converge a cero 14 2.3.1 Correspondencia con el problema de optimización con matrices no diagonales 15 2.4 Control Predictivo Basado en Modelo que converge a una referencia 17 2.4.1 Correspondencia con el problema de optimización con restricciones de igualdad 18 2.5 Steady-State Target Optimization 21 2.5.1 Filtro del SSTO 23 2.5.2 Ejemplo de aplicación sobre la planta de los 4 tanques 23 2.6 Observador de estado 25 2.6.1 Filtro del observador de estado 26 2.6.2 Ejemplo de aplicación sobre la planta de los 4 tanques 26 3 Implementación en LabVIEW 29 3.1 Programación del Control Predictivo Basado en Modelo 29 3.2 Aspectos a tener en cuenta antes de la implementación 30 3.2.1 Simulación secuencial previa 30 3.2.2 Muestreo en paralelo 32 3.2.3 Control Predictivo de la planta de los 4 tanques 34 Referencias y Bibliografía 37 Glosario 39 Anexo A: Scripts de MATLAB 41 Anexo B: DLL 49 xvii ÍNDICE DE TABLAS Tabla 1-1. Solución del problema 1 2 Tabla 1-2. Solución del problema 2 3 Tabla 2-1. Puntos de funcionamiento de la planta de los 4 tanques 24 xix ÍNDICE DE FIGURAS Figura 1-1. Parábola convexa 2 Figura 1-2. Solución del problema 1 2 Figura 1-3. Comparativa de tiempos de convergencia QUADPROG – ISTA 6 Figura 1-4. Comparativa de tiempos de convergencia QUADPROG – FISTA 7 Figura 1-5. Comparativa de tiempos de convergencia con restricciones de igualdad 10 Figura 2-1. Control Predictivo Basado en Modelo que converge a cero 17 Figura 2-2. Control Predictivo Basado en Modelo que converge a una referencia 20 Figura 2-3. Control Predictivo Basado en Modelo que converge a una referencia para 𝑝 = 𝑚 21 Figura 2-4. Diagrama de bloques SSTO-MPC 23 Figura 2-5. Planta de Johansson 24 Figura 2-6. Simulación con el modelo lineal de la planta de los 4 tanques 25 Figura 2-7. Diagrama de bloques SSTO-MPC-Observador 26 Figura 2-8. Simulación con el modelo no lineal de la planta de los 4 tanques 27 Figura 2-9. Estado real y observado para la simulación con el modelo no lineal 27 Figura 2-10. Diagrama de bloques MPC 28 Figura 3-1. Block Diagram de la simulación secuencial previa en LabVIEW 31 Figura 3-2. Comparativa de las simulaciones secuenciales MATLAB – LabVIEW 31 Figura 3-3. Diferencia entre las simulaciones secuenciales MATLAB – LabVIEW 32 Figura 3-4. Muestreo ideal de la planta de los 4 tanques 33 Figura 3-5. Muestreo real de la planta de los 4 tanques 33 Figura 3-6. Block Diagram del Control Predictivo Basado en Modelo 34 Figura 3-7. Block Diagram del modelo no lineal de la planta de los 4 tanques 35 Figura 3-8. Front Panel 36 Figura 3-9. Simulación MPC – Modelo no lineal de la planta de los 4 tanques 36 xxi Notación : Tal que < Menor > Mayor ≥ Menor o igual ≤ Mayor o igual s.*a.* Sujeto a ∈* Perteneciente ∆ Incremento ∇ Jacobiano 𝐴0 Transpuesta 𝐼 Matriz identidad 𝑧 Valor absoluto 𝑧∗* Óptimo 𝑥5 Valor inicial 𝑥6 Valor en el instante de tiempo 𝑘 𝑥 Valor estimado 1 1 INTRODUCCIÓN A LA OPTIMIZACIÓN l objetivo de la optimización consiste en calcular el conjunto de variables que minimizan (o maximizan) un determinado criterio satisfaciendo un conjunto de restricciones sobre las mismas. El criterio se debe elegir de manera adecuada para hacer que la solución sea única. En este capítulo nos vamos a centrar en resolver problemas de optimización cuyo coste es cuadrático. Para la resolución de estos problemas van a cobrar gran importancia los algoritmos que se utilicen para resolverlos. Éstos deben ser sencillos de implementar en computadores de bajas prestaciones para así poder ser utilizados por un amplio abanico de equipos de control. 1.1 Problema 1 Sea 𝑧 un número real, calcular 𝑧∗, la solución del problema min ;**********1 2𝑑𝑧?+𝑐𝑧 [1.1] 𝑠. 𝑎.***********𝑎 ≤ 𝑧 ≤ 𝑏 [1.2] siendo 𝑎 < 𝑏 y 𝑑 > 0. Solución La función para la que buscamos la 𝑧 mínima, se corresponde con la ecuación de una parábola definida tal que 𝑦 = 1 2𝑑 𝑧 − 𝑧5?+ 𝑦5 [1.3] E Introducción a la optimización 2 Figura 1-1. Parábola convexa. Es fácil comprobar gráficamente que la 𝑧 que hace mínima 𝑦 será 𝑧5 siempre que ésta esté comprendida entre los valores 𝑎 y 𝑏. En caso contrario, tomará el valor del límite más cercano tal que Tabla 1-1. Solución del problema 1. Si 𝑎 ≤ 𝑧5≤ 𝑏 𝑧∗= 𝑧5 Si 𝑧5> 𝑏 𝑧∗= 𝑏 Si 𝑧5< 𝑎 𝑧∗= 𝑎 Figura 1-2. Solución del problema 1. 3 3 Implementación de controladores predictivos en LabVIEW Como estamos buscando la solución mínina, al desarrollar la ecuación [1.3], se pueden eliminar todos los términos constantes que no sean dependientes de 𝑧. 𝑦 = 1 2𝑑𝑧?+1 2𝑑𝑧5 ?−𝑑𝑧𝑧5+ 𝑦5 [1.4] 𝑦 = 1 2𝑑𝑧?−𝑑𝑧𝑧5 [1.5] Por lo tanto, igualando [1.5] con la función [1.1] que queremos minimizar, obtenemos 1 2𝑑𝑧?−𝑑𝑧𝑧5=1 2𝑑𝑧?+𝑐𝑧 [1.6] 𝑧5= − 𝑐 𝑑 [1.7] 1.2 Problema 2 Sea 𝑧 ∈ ℝI un vector de números reales, calcular 𝑧∗, la solución del problema min ;**********1 2𝑧0𝐷𝑧 + 𝐶0𝑧 [1.8] 𝑠. 𝑎.***********𝑎L≤ 𝑧L≤ 𝑏L,*****𝑖 = 1, … , 𝑛 [1.9] siendo 𝑎L< 𝑏L y siendo 𝐷 ∈ ℝI*×*I una matriz diagonal tal que 𝑑LL > 0. Solución Se comprueba que la solución de este problema es similar a la del problema anterior. Desarrollando la ecuación [1.8] se obtiene 1 2𝑧R ?𝑑RR + 𝑐R𝑧R+ ⋯ + 1 2𝑧I ?𝑑II + 𝑐I𝑧I [1.10] Por lo que nuestro problema se puede separar en 𝑛 problemas distintos, siendo cada uno de ellos, similar al problema 1. La solución será 𝑧∗∈ ℝI, cuyas componentes se corresponden con las soluciones individuales de cada problema. Tabla 1-2. Solución del problema 2. Si 𝑎L≤ 𝑧5T≤ 𝑏L 𝑧L ∗= 𝑧5T Si 𝑧5T> 𝑏L 𝑧L ∗= 𝑏L Si 𝑧5T< 𝑎L 𝑧L ∗= 𝑎L Donde 𝑧5T*= − 𝑐L 𝑑LL [1.11] Introducción a la optimización 10 convergencia del algoritmo FISTA con la función quadprog. Una vez más, observamos que el algoritmo FISTA presenta una clara ventaja en tiempo de convergencia. Figura 1-5. Comparativa de tiempos de convergencia con restricciones de igualdad. 11 2 INTRODUCCIÓN AL CONTROL PREDICTIVO l concepto de Control Predictivo Basado en Modelo, en inglés, Model Predictive Control, se basa en el uso de un modelo dinámico para predecir el comportamiento del sistema y optimizar dicha predicción para poder realizar el control en tiempo real. Cobran gran importancia los modelos del sistema ya que el control óptimo del sistema depende del estado inicial del mismo. 2.1 Modelos de sistemas El modelo dinámico de un sistema queda definido por las siguientes ecuaciones diferenciales 𝑑𝑋 𝑑𝑡 = 𝐹(𝑋, 𝑈) [2.1] 𝑌 = 𝐻(𝑋, 𝑈) [2.2] 𝑋 0 = 𝑋5 [2.3] donde 𝑋 ∈ ℝI representa el vector de estados del sistema, 𝑌 ∈ ℝZ el vector de entradas modificables del sistema, 𝑌 ∈ ℝ‚ el vector de salidas controlables del sistema y 𝑋5∈ ℝI el estado inicial del sistema. 2.1.1 Modelo dinámico lineal invariante con el tiempo Este modelo es linealizable y además, se puede representar sin que varíe con el tiempo. Si linealizamos el sistema en un punto de funcionamiento (𝑋ƒ, 𝑈ƒ) el modelo dinámico lineal invariante con el tiempo queda expresado como 𝑑𝑥 𝑑𝑡 =𝐴𝑥 +𝐵𝑢 [2.4] 𝑦 = 𝐶𝑥 +𝐷𝑢 [2.5] 𝑥 0 = 𝑥5 [2.6] E Introducción al Control Predictivo 12 siendo 𝐴 ∈ ℝI×I 𝐴 = ∇\𝐹(𝑋, 𝑈) (…†,y†) [2.7] la matriz de transición de estados y 𝐵 ∈ ℝI×Z 𝐵 = ∇‡𝐹(𝑋, 𝑈) (…†,y†) [2.8] la matriz de entrada siendo 𝐴 𝐵 el jacobiano de 𝐹. La matriz 𝐶 ∈ ℝ‚×I representa la matriz de salida y 𝐷 ∈ ℝ‚×Z la matriz que relaciona directamente la entrada modificable con la salida controlable. El estado del sistema 𝑥 ∈ ℝI, la entrada modificable 𝑢 ∈ ℝZ y la salida controlable 𝑦 ∈ ℝZ están expresados en variables incrementales con respecto al punto de funcionamiento sobre el que hemos linealizado de manera que 𝑥 = 𝑋 − 𝑋ƒ [2.9] 𝑢 = 𝑈 − 𝑈ƒ [2.10] 𝑦 = 𝑌 − 𝑋ƒ [2.11] siendo 𝑥5∈ ℝI el estado inicial del sistema expresado en variables incrementales. 2.1.2 Modelo dinámico lineal invariante con el tiempo discretizado Este modelo dinámico lineal invariante con el tiempo se puede expresar discretizado en el tiempo. Si se muestrea para un cierto periodo T, el modelo dinámico lineal invariante en tiempo discreto se expresa como 𝑥6dR = 𝐴𝑥6+ 𝐵𝑢6 [2.12] 𝑦6= 𝐶𝑥6+ 𝐷𝑢6 [2.13] 𝑥5= 𝑥5 [2.14] 2.1.3 Restricciones en tiempo Las entradas y salidas de los sistemas suelen estar limitadas. Es necesario por tanto, limitar el estado 𝑥 y la actuación 𝑢 del sistema. Estas restricciones se denominan restricciones en tiempo y se deben cumplir para todo 𝑘, quedando éstas expresadas como 𝑥6∈ 𝑋 [2.15] 𝑢6∈ 𝑈 [2.16] siendo 𝑋 = 𝑥:*****𝐿𝐵\≤ 𝑥 ≤ 𝑈𝐵\ [2.17] 𝑈 = 𝑢:*****𝐿𝐵‡≤𝑢≤𝑈𝐵‡ [2.18] donde los vectores 𝐿𝐵\ y 𝑈𝐵\ de dimensión ℝI representan los límites inferiores (Lower Bounds) y superiores (Upper Bounds) de los 𝑛 estados del sistema, y los vectores 𝐿𝐵‡ y 𝑈𝐵‡ de dimensión ℝZ representan los límites inferiores y superiores de las 𝑚 entradas manipulables del sistema. 13 13 Implementación de controladores predictivos en LabVIEW 2.2 El controlador óptimo El objetivo del controlador óptimo es que la salida controlable del sistema se estabilice en torno a una referencia 𝑟 ∈ ℝ‚ tal que 𝑦6= 𝑟 [2.19] Para ello, el sistema se debe estabilizarse en un punto de equilibrio en el que se cumpla 𝑥‰= 𝐴𝑥‰+ 𝐵𝑢‰ [2.20] 𝑟 = 𝐶𝑥‰+ 𝐷𝑢‰ [2.21] Una vez definidos los siguientes cambios de variables 𝑥6= 𝑥6− 𝑥‰ [2.22] 𝑢6= 𝑢6− 𝑢‰ [2.23] el objetivo del controlador óptimo pasa a ser 𝑥 → 0 [2.24] 𝑢 → 0 [2.25] La función de coste cuadrático es definida como 𝐿 𝑥, 𝑢 = 𝑥0𝑄𝑥 + 𝑢0𝑅𝑢 [2.26] donde 𝑄 ∈ ℝI y 𝑅 ∈ ℝZ. Esta función mide el rendimiento del sistema a partir de la entrada manipulable 𝑢 y el estado 𝑥 del sistema. El coste de la trayectoria predicha para un horizonte de predicción 𝑁 se puede expresar como 𝑉•𝑥, 𝑢 = 𝐿(𝑥6, 𝑢6) • 6Ž5 [2.27] Entonces, la ley de control óptima queda definida por min ‡*****𝑉 •𝑥, 𝑢 = 𝑥6 0𝑄𝑥6+ 𝑢6 0𝑅𝑢6 • 6Ž5 [2.28] 𝑥5= 𝑥5 [2.29] 𝑥6dR = 𝐴𝑥6+ 𝐵𝑢6 [2.30] 𝐿𝐵\≤ 𝑥6≤𝑈𝐵\ [2.31] 𝐿𝐵‡≤ 𝑢6≤𝑈𝐵‡ [2.32] siendo 𝐿𝐵\=𝐿𝐵\− 𝑥‰ [2.33] 𝑈𝐵\=𝑈𝐵\− 𝑥‰ [2.34] Introducción al Control Predictivo 14 𝐿𝐵‡=𝐿𝐵‡− 𝑢‰ [2.35] 𝑈𝐵‡=𝑈𝐵‡− 𝑢‰ [2.36] La ley de control óptima definida calcula el coste de la trayectoria del sistema para un horizonte de predicción infinito. El controlador óptimo tiene como función minimizar esta trayectoria con horizonte infinito. 2.2.1 Ley de control del MPC El controlador óptimo no tiene una solución explícita y en la práctica solo puede resolverse para ciertos casos como para sistemas sin restricciones. La ley de control del Control Predictivo Basado en Modelo se basa en la ley de control óptima para un horizonte de predicción 𝑁 < ∞. El cálculo de la trayectoria del sistema para un horizonte de predicción 𝑁 definido tiene asociada una restricción adicional llamada restricción terminal. Esta restricción implica que el estado terminal del sistema 𝑥• tiene que estar comprendido dentro de unos límites 𝑋‘∈ ℝI. Por tanto, la ley de control del Control Predictivo Basado en Modelo para un horizonte de predicción 𝑁 finito quedaría expresada como min ‡𝑥6 0𝑄𝑥6+ 𝑢6 0𝑅𝑢6 •mR 6Ž5 + 𝑥• 0𝑃𝑥• [2.37] 𝑥5= 𝑥5 [2.38] 𝑥6dR = 𝐴𝑥6+ 𝐵𝑢6 [2.39] 𝐿𝐵\≤ 𝑥6≤𝑈𝐵\ [2.40] 𝐿𝐵‡≤ 𝑢6≤𝑈𝐵‡ [2.41] 𝑥•∈ 𝑋‘ [2.42] siendo 𝑃 ∈ ℝI y la solución del problema 𝑢∗= 𝑢5 ∗ ⋮ 𝑢•mR ∗ [2.43] 2.2.1.1 Ley de control De la solución óptima del problema de optimización se deriva la ley de control gracias a la estrategia de horizonte deslizante. Así se aplica sólo la acción de control óptima actual, recalculándose el problema de optimización en el instante siguiente para el nuevo estado del sistema. El MPC solo aplica la actuación 𝑢5 ∗ calculada en el instante de tiempo 𝑘 = 0. Para el siguente instante de tiempo calcula de nuevo la solución del problema usando el resto de la solución actual 𝑢‘ ∗ como solución inicial del problema. De nuevo, aplica únicamente 𝑢5 ∗ y recalcula la solución del problema así sucesivamente. Al realizar el control en bucle cerrado los resultados son mejores ya que, al realimentar el controlador predictivo, obtiene el estado real del sistema para el siguiente instante de tiempo 𝑘. 2.3 Control Predictivo Basado en Modelo que converge a cero La ley de control del Control Predictivo Basado en Modelo para un sistema lineal discreto invariante con el tiempo que converge a cero de manera que 𝑟 = 0 se puede expresar como 15 15 Implementación de controladores predictivos en LabVIEW min ‡𝑥6 0𝑄𝑥6+ 𝑢6 0𝑅𝑢6 •mR 6Ž5 [2.44] 𝑥5= 𝑥5 [2.45] 𝑥6dR = 𝐴𝑥6+ 𝐵𝑢6 [2.46] 𝐿𝐵\≤ 𝑥6≤𝑈𝐵\ [2.47] 𝐿𝐵‡≤ 𝑢6≤𝑈𝐵‡ [2.48] 𝑥•= 0 [2.49] 2.3.1 Correspondencia con el problema de optimización con matrices no diagonales Este problema se puede expresar como el problema de optimización con restricciones linealmente independientes resuelto en el capítulo anterior, quedando éste expresado como min ;**********𝑧0𝐻𝑧 [2.50] 𝐸𝑧 = 𝑏 [2.51] 𝐿𝐵;≤ 𝑧 ≤ 𝑈𝐵; [2.52] Siendo 𝑧 ∈ ℝ•∗(IdZ), de la forma 𝑧 = 𝑥5 𝑢5 ⋮ 𝑥•mR 𝑢•mR [2.53] Siendo 𝐻 ∈ ℝ•∗ IdZ *×*•∗(IdZ), de la forma 𝐻 = 𝑄 0 0 𝑅 ⋯0 0 0 0 ⋮ ⋱ ⋮ 0 0 0 0 ⋯𝑄 0 0 𝑅 [2.54] Siendo 𝐸 ∈ ℝI∗ •dR *×*•∗(IdZ), de la forma 𝐸 = 𝐼 0 𝐴 𝐵 0 0 000 −𝐼 0 0 𝐴 𝐵 −𝐼 ⋯ 0 0 0 0 0 0 ⋮ ⋱ ⋮ 00000 ⋯𝐴 𝐵 [2.55] Siendo 𝑏 ∈ ℝI∗(•dR), de la forma 𝑏 = 𝑥5 0 ⋮ 0 [2.56] Introducción al Control Predictivo 16 Siendo 𝐿𝐵;∈ ℝ•∗(IdZ), de la forma 𝐿𝐵;= 𝐿𝐵\ 𝐿𝐵‡ ⋮ 𝐿𝐵\ 𝐿𝐵‡ [2.57] Siendo 𝑈𝐵;∈ ℝ•∗(IdZ), de la forma 𝑈𝐵;= 𝑈𝐵\ 𝑈𝐵‡ ⋮ 𝑈𝐵\ 𝑈𝐵‡ [2.58] Ejemplo 1 A modo de comprobación, creamos con MATLAB un sistema de dimensión ℝ”×? constituido por 3 estados, 2 entradas modificables y 2 salidas controlables. Esto implica la creación de las matrices 𝐴, 𝐵, 𝐶 y 𝐷, las matrices 𝑄 y 𝑅, los límites 𝐿𝐵\, 𝐿𝐵‡, 𝑈𝐵\ y 𝑈𝐵‡ así como el vector de estados iniciales 𝑥5, siendo estos 𝐴 = 0.8 0.1 0 0.1 0.7 0 0 0 0.85 𝐵 = 0.5 0.1 1 0 0.1 0.7 𝐶 = 1 0 0 0 1 1 𝐷 = 0 0 0 0 𝑄 = 50 0 0 050 0 0 0 50 𝑅 = 10 0 010 𝐿𝐵\= −2 −2 −2 𝐿𝐵‡=−3 −3 𝑈𝐵\= 2 2 2 𝑈𝐵‡=3 3 17 17 Implementación de controladores predictivos en LabVIEW 𝑥5= 1.68 1.56 1.6 A continuación se aplica en bucle cerrado el Control Predictivo Basado en Modelo resolviéndolo mediante el algoritmo de optimización FISTA para un horizonte de predicción 𝑁 = 10 y durante 20 intervalos de muestreo. En cada intervalo de muestreo se toma, de la solución 𝑧∗ del problema, el estado 𝑥R que resultaría de aplicar la actuación 𝑢5 al sistema estando el sistema en el estado 𝑥5. Este estado 𝑥R será el nuevo estado 𝑥5 en el siguiente intervalo de muestreo. Al representar el vector 𝑥5 para cada intervalo de nuestreo se observa la convergencia de los estados del sistema a cero. Figura 2-1. Control Predictivo Basado en Modelo que converge a cero. 2.4 Control Predictivo Basado en Modelo que converge a una referencia Para un sistema lineal discreto invariante con el tiempo, la ley de control óptima para una convergencia en torno a una referencia 𝑟 se puede expresar como min ‡𝑥6− 𝑥‰0𝑄(𝑥6− 𝑥‰) + 𝑢6− 𝑢‰0𝑅(𝑢6− 𝑢‰) •mR 6Ž5 [2.59] 𝑥5= 𝑥5 [2.60] 𝑥6dR = 𝐴𝑥6+ 𝐵𝑢6 [2.61] 𝐿𝐵\≤ 𝑥6≤𝑈𝐵\ [2.62] 𝐿𝐵‡≤ 𝑢6≤𝑈𝐵‡ [2.63] 𝑥•= 𝑥‰ [2.64] Introducción al Control Predictivo 18 2.4.1 Correspondencia con el problema de optimización con restricciones de igualdad De nuevo, este problema de minimización puede volver a expresarse como el problema de optimización con restricciones de igualdad desarrollado en el capítulo anterior. min ;**********𝑧0𝐻𝑧 + 𝑓0𝑧 [2.65] 𝐸𝑧 = 𝑏 [2.66] 𝐿𝐵;≤ 𝑧 ≤ 𝑈𝐵; [2.67] Siendo 𝑧 ∈ ℝ•∗(IdZ) de nuevo de la forma 𝑧 = 𝑥5 𝑢5 ⋮ 𝑥•mR 𝑢•mR [2.68] Siendo 𝐻 ∈ ℝ•∗ IdZ *×*•∗(IdZ), de la forma 𝐻 = 𝑄 0 0 𝑅 ⋯0 0 0 0 ⋮ ⋱ ⋮ 0 0 0 0 ⋯𝑄 0 0 𝑅 [2.69] Siendo 𝑓 ∈ ℝ•∗(IdZ), de la forma 𝑓 = −𝑄𝑥‰ −𝑅𝑢‰ ⋮ −𝑄𝑥‰ −𝑅𝑢‰ [2.70] Siendo 𝐸 ∈ ℝI∗ •dR *×*•∗(IdZ), de la forma 𝐸 = 𝐼 0 𝐴 𝐵 0 0 000 −𝐼 0 0 𝐴 𝐵 −𝐼 ⋯ 0 0 0 0 0 0 ⋮ ⋱ ⋮ 00000 ⋯𝐴 𝐵 [2.71] Siendo 𝑏 ∈ ℝI∗(•dR), de la forma 𝑏 = 𝑥5 0 ⋮ 𝑥‰ [2.72] Siendo 𝐿𝐵;∈ ℝ•∗(IdZ), de la forma 𝐿𝐵;= 𝐿𝐵\ 𝐿𝐵‡ ⋮ 𝐿𝐵\ 𝐿𝐵‡ [2.73] 19 19 Implementación de controladores predictivos en LabVIEW Siendo 𝑈𝐵;∈ ℝ•∗(IdZ), de la forma 𝑈𝐵;= 𝑈𝐵\ 𝑈𝐵‡ ⋮ 𝑈𝐵\ 𝑈𝐵‡ [2.74] Ejemplo 2 Aplicamos de nuevo sobre un sistema de dimensión ℝ”×? creado con MATLAB el Control Predictivo Basado en Modelo en bucle cerrado resuelto mediante el algoritmo de optimización FISTA. En este caso, a las matrices 𝐴, 𝐵, 𝐶 y 𝐷, las matrices 𝑄 y 𝑅, los límites 𝐿𝐵\, 𝐿𝐵‡, 𝑈𝐵\ y 𝑈𝐵‡ y el vector de estados iniciales 𝑥5 deben añadirse ahora, los vectores 𝑥‰ y 𝑢‰. Las ecuaciones [2.20] y [2.21] pueden ser expresadas tal que 𝐴 − 𝐼 𝐵 𝐶 𝐷 𝑥‰ 𝑢‰=0 𝑟 [2.75] De manera que el cálculo de los vectores 𝑥‰ y 𝑢‰ se realizaría tal que 𝑥‰ 𝑢‰=𝐴 − 𝐼 𝐵 𝐶 𝐷 \0 𝑟 [2.76] Por lo tanto, las matrices y vectores genearados son 𝐴 = 0.8 0.1 0 0.1 0.7 0 0 0 0.85 𝐵 = 0.5 0.1 1 0 0.1 0.7 𝐶 = 100 011 𝐷 = 0 0 0 0 𝑄 = 50 0 0 050 0 0 0 50 𝑅 = 10 0 010 𝐿𝐵\= −20 −20 −20 𝐿𝐵‡=−30 −30 𝑈𝐵\= 20 20 20 Introducción al Control Predictivo 26 y calculado tal que 𝐿5:* 𝜆L(𝐴5− 𝐿5𝐶5) < 1 [2.104] Siendo 𝐶5∈ ℝZ*×*(IdZ), de la forma 𝐶5=𝐶 𝐼 [2.105] 2.6.1 Filtro del observador de estado Al igual que el SSTO, el observador tiene asociado un filtro con el fin de evitar cambios bruscos al producirse cambios en la perturbación del sistema. Por lo tanto, la perturbación medida por el observador estará definidapor una 𝑑‰ calculada a partir de 𝑑‰l~• = 1 − 𝛼 𝑑‰l+ 𝛼𝑑 [2.106] siendo 𝛼 ∈ [0, 1] y 𝑑‰ la perturbación estimada. 2.6.2 Ejemplo de aplicación sobre la planta de los 4 tanques Realizamos en MATLAB una simulación para comprobar el buen funcionamiento del MPC, el SSTO y el observador de estado conjuntamente según el siguiente diagrama de bloques. Figura 2-7. Diagrama de bloques SSTO-MPC-Observador La perturbación estimada a la salida del sistema 𝑑‰ debe ser realimentada hacia la referencia de manera que 𝑟 = 𝑟 − 𝑑‰ [2.107] Para para poder generar una perturbación a la salida nuestro sistema va a ser el modelo no lineal de la planta de los 4 tanques de Johansson. El modelo de la planta que tiene el MPC es el linealizado en torno a un estado 𝑋ƒ. Esto hace que exista una diferencia entre el sistema que vamos a controlar y el modelo contenido en el MPC cuando el estado del sistema sea distinto al del punto de linealización 𝑋ƒ. Esta diferencia será interpretada por el observador como 𝑑. En la simulación realizada en MATLAB haremos un cambio de referencia en el intervalo de muestreo número 500. De igual manera, cada 200 intervalos de muestreo, añadiremos una perturbación aleatoria a la salida controlable del sistema. 27 27 Implementación de controladores predictivos en LabVIEW Figura 2-8. Simulación con el modelo no lineal de la planta de los 4 tanques Se puede observar cómo el sistema converge a la referencia con exactitud a pesar de que el MPC está diseñado para un sistema distinto al que está controlando. Además, es capaz de llevar la salida controlable a la referencia aunque ésta tenga una perturbación adicional. Es interesante observar cómo el estado del sistema 𝑋 y el estado del sistema estimado por el observador 𝑋 no coinciden. Esto se debe a la perturbación que afecta a la salida controlable 𝑌 ya sea por una perturbación o por tratarse de un modelo no lineal. El observador estima la perturbación que afecta a la salida 𝑑 y el estado en el que tendría que estar el sistema para contrarrestar esa perturbación estimada y encontrarse en la referencia deseada. La perturbación estimada a la salida del sistema es también restada a la referencia de entrada para conseguir llevar al sistema a un estado distinto en el que junto a la perturbación, alcance la referencia correctamente. Figura 2-9. Estado real y observado para la simulación con el modelo no lineal El diagrama de bloques de la simulación puede ser también expresado de una manera más compacta incluyendo dentro del bloque del MPC el observador, el SSTO y sus correspondientes filtros. Además se Introducción al Control Predictivo 28 incluyen los cambios de variables absolutas a variables incrementales y viceversa para el correcto funcionamiento conjunto del MPC y el sistema. Figura 2-10. Diagrama de bloques MPC 29 3 IMPLEMENTACIÓN EN LABVIEW abVIEW es una plataforma para el diseño de sistemas especializada en adquisición de datos, instrumentación y control, y automatización industrial. El lenguaje de programación usado por LabVIEW es visual y se basa en el fujo de datos. La ejecución del mismo se determina a través de la estructura de los diagramas de bloques. 3.1 Programación del Control Predictivo Basado en Modelo Para implementar el Control Predictivo Basado en Modelo en LabVIEW utilizaremos una DLL. Se trata de una biblioteca de enlace dinámico que tiene como función almacenar código que puede ser ejecutado. En esta biblioteca estará incluida nuestra función MPC, programada en lenguaje C, que contiene al controlador predictivo. Por comodidad, a través de un script de MATLAB se genera un archivo de texto .txt en el que se definen todas las variables y constantes necesarias para la cabecera de la función MPC. De manera que al cambiar los parámetros modificables del Control Predictivo Basado en Modelo en el script de MATLAB, se genere automáticamente la cabecera correcta de la función para los parámetros indicados. Además dicho script se encargará también de realizar cálculos previos con los parámetros modificables del controlador de manera que la función MPC solo tenga que realizar aquellos cálculos que dependan exclusivamente de las entradas de dicha función. De esta manera se evita que la función repita ciertos cálculos que siempre dan el mismo resultado cada vez que se ejecute dicha función. Del mismo modo, para hacer más rápida la ejecución de la función, algunas matrices y vectores que sigan un patrón de repetición o contengan una gran cantidad de ceros será definidas para que tengan una dimensión menor. Las matrices 𝑄 y 𝑅 serán redefinidas como los vectores 𝑄 ∈ ℝI y 𝑅 ∈ ℝZ. La matriz 𝐸 pasa a estar definida como 𝐸 ∈ ℝI*×*IdZ tal que 𝐸 = 𝐴 𝐵 −𝐼 [3.1] Para usar la matriz 𝐸 de manera traspuesta se definen dos matrices adicionales 𝐸𝑇1 ∈ ℝIdZ*×*?∗I tal que 𝐸𝑇1 = 𝐼 𝐴 0 𝐵 [3.2] L Implementación en LabVIEW 30 y 𝐸𝑇2 ∈ ℝIdZ*×*?∗I definida como 𝐸𝑇2 = −𝐼 𝐴 0 𝐵 [3.3] El vector 𝑓 basta con estar definido como 𝑓 ∈ ℝIdZ 𝑓 = −𝑄𝑥‰ −𝑅𝑢‰ [3.4] Los vectores 𝑈𝐵\, 𝐿𝐵\, 𝑈𝐵‡ y 𝐿𝐵‡ pasan a estar definidos por los vectores 𝑒𝑑𝑔𝑒 y 𝑛𝑒𝑤𝑧 que representan un cálculo intermedo del Control Predictivo Basado en Modelo que se repite para cada ejecución del mismo. Siendo 𝑒𝑑𝑔𝑒 ∈ ℝIdZ definido tal que 𝑒𝑑𝑔𝑒 = 𝑈𝐵\−𝐿𝐵\ 2 𝑈𝐵‡−𝐿𝐵‡ 2 [3.5] y 𝑛𝑒𝑤𝑧 ∈ ℝIdZ como 𝑛𝑒𝑤𝑧 = 𝑈𝐵\+𝐿𝐵\ 2 𝑈𝐵‡+𝐿𝐵‡ 2 [3.6] 3.2 Aspectos a tener en cuenta antes de la implementación La principal diferencia entre las simulaciones que se van a realizar en LabVIEW y las programadas hasta ahora en MATLAB es que éstas son secuenciales mientras que en LabVIEW las tareas de control y simulación o muestreo de la planta se ejecutan en paralelo de acuerdo a la precedencia impuesta por LabVIEW. Esto hace que no se puedan comparar directamente dos simulaciones iguales realizadas en MATLAB y en LabVIEW. Además, es interesante hacer un estudio previo sobre el muestreo en tiempo real que se va a realizar en LabVIEW. 3.2.1 Simulación secuencial previa En primer lugar se va a comparar la correcta implementación en lenguaje C de la función MPC. Para ello, se realiza una simulación secuencial en LabVIEW de manera que esta simulación sea similar a una simulación programada en MATLAB. Dentro de una estructura While Loop que se repite cada 500 milisegundos se añade una estructura Flat Sequence en la que se simularán de forma secuencial la planta de los 4 tanques y el controlador predictivo. En el primer subdiagrama de esta estructura se incluye un bloque Call Library Function que cargará la función 𝑃𝐿𝐴𝑁𝑇𝐴(𝑈, 𝑋, 𝑌) contenida en una DLL. Esta función simulará el modelo no lineal de la planta de los 4 tanques para un tiempo de muestreo. A la salida de la misma, se tomarán como datos la salida controlable 𝑌 del sistema. En el segundo subdiagrama de la estructura Flat Sequence se incluye otro bloque Call Library Function que llama a la función 𝑀𝑃𝐶(𝑈, 𝑌, 𝑅𝑒𝑓) contenida en la misma DLL anterior. Para esta simulación previa, el estado estimado por el observador 𝑋, la perturbación estimada por el observador 𝑑, la referencia filtrada 𝑟 ƒL¡k‰[¢[ y la perturbación estimada filtrada 𝑑‰ serán definidos como stactic dentro de la función y no serán accesibles. 31 31 Implementación de controladores predictivos en LabVIEW Figura 3-1. Block Diagram de la simulación secuencial previa en LabVIEW Como se puede observar a simple vista, los resultados de ambas simulaciones coindicen notablemente para ambas salidas. Figura 3-2. Comparativa de las simulaciones secuenciales MATLAB – LabVIEW Para una mejor comparativa, se muestra la diferencia en valor absoluto entre los resultados obtenidos de la simulación programada en MATLAB y los resultados obtenidos de la simulación secuencial en LabVIEW. Implementación en LabVIEW 32 Figura 3-3. Diferencia entre las simulaciones secuenciales MATLAB – LabVIEW La pequeña diferencia que se observa se debe a la precisión con la que realizan los cálculos MATLAB y LabVIEW. Mientras que en MATLAB las variables necesarias para el MPC son deducidas unas a partir de otras y acarreando todos los decimales, en LabVIEW todas las variables han sido definidas en la cabecera de la función con un máximo de 4 cifras decimales. 3.2.2 Muestreo en paralelo Antes de implementar por completo el Control Predictivo Basado en Modelo en LabVIEW es conveniente hacer un estudio sobre el muestreo de la planta de los 4 tanques en paralelo. Para ello, se simula únicamente el modelo no lineal de la planta de los 4 tanques de la siguiente manera: - Cada 100 ms se ejecuta la función 𝑃𝐿𝐴𝑁𝑇𝐴(𝑈, 𝑋, 𝑌) que ejecuta el modelo no lineal de la planta una vez y se toma como dato la lectura de la salida controlable 𝑌. - Cada 300 ms se incrementa la entrada modificable 𝑈 en 0.01. - Cada 5000 ms se toma como dato muestreado la lectura de la salida controlable 𝑌. Para esta simulación programada en MATLAB y en LabVIEW, el muestreo se realiza idealmente en los instantes de tiempo 151 101 151 201 ⋯ Es decir, se muestrea la planta de los 4 tanques una vez cada 50 veces que se ejecuta la función que contiene su modelo no lineal. 33 33 Implementación de controladores predictivos en LabVIEW Figura 3-4. Muestreo ideal de la planta de los 4 tanques Si observamos el muestreo realizado en MATLAB, coincide con la lectura de la planta para esos instantes de tiempo. Sin embargo, para el muestreo realizado en LabVIEW, los valores muestreados no coindicen con los instantes de tiempo indicados. Figura 3-5. Muestreo real de la planta de los 4 tanques Tras un análisis exhaustivo sobre los datos muestreados en LabVIEW, se observa que los valores muestreados de la planta coinciden con la planta simulada en los instantes de tiempo 144 94 144 194 ⋯ Se deduce por tanto que en paralelo no se realiza correctamente el muestreo de la planta durante los primeros instantes de tiempo, pero sin embargo, tras varios muestreos incorrectos, sí realiza correctamente el muestreo cada 50 veces que se ejecuta el modelo no lineal de la planta. Como solución a este problema, se propone hacer que la simulación de la planta de los 4 tanques parta siempre Implementación en LabVIEW 34 de un punto de funcionamiento manteniendo la referencia en dicho punto. En ese caso la toma de datos en los primeros instantes de muestreo se hará de forma correcta ya que mientras se mantenga la referencia, el valor de la salida no varía al tratarse de un punto de funcionamiento de la planta. 3.3 Control Predictivo de la planta de los 4 tanques Para la simulación completa del Control Predictivo Basado en Modelo en LabVIEW, se implementará el controlador dentro de una estructura While Loop. Este bucle se ejecutará cada 5000 milisegundos siendo este el tiempo de muestreo del controlador. Para hacer que la simulación sea en tiempo real para el tiempo de muestreo indicado se incluye la estructura Wait Until Next ms Multiple dentro del bucle while. Esto permitirá que se ejecute el bucle en tiempo real para cada tiempo de muestreo independiente de la duración de la ejecución del contenido del bucle. La función MPC contenida en el archivo .dll será llamada durante la simulación mediante el bloque Call Library Function Node. En este bloque es necesario configurar la ruta donde se encuentra la DLL, el nombre de la función y el prototipo de la función, siendo ésta 𝑀𝑃𝐶(𝑈, 𝑌, 𝑅𝑒𝑓,𝑂𝑏𝑠,𝑟𝑓𝑖𝑙𝑡𝑟𝑎𝑑𝑎,𝑑𝑔𝑜𝑟𝑟𝑜,𝑑𝑟) Siendo el argumento 𝑈 la actuación calculada por el controlador, el argumento 𝑌 la salida muestreada de la planta, 𝑅𝑒𝑓 la referencia alcanzable, 𝑂𝑏𝑠 el estado estimado por el observador, 𝑟𝑓𝑖𝑙𝑡𝑟𝑎𝑑𝑎 la referencia filtrada, 𝑑𝑔𝑜𝑟𝑟𝑜 la perturbación estimada por el observador y 𝑑𝑟 la perturbación estimada filtrada. Estos últimos 4 argumentos deben ser inicializados junto a la entrada 𝑈 inicial. El tipo de datos de todos ellos serán array data pointer con una resolución double de 8 byte. La entrada modificable 𝑈 y salida controlable 𝑌 de la planta utilizadas en la función MPC se definirán mediante estructuras Global Variable. Esta estructura hace que dichas variables sean de tipo global. La referencia por su parte, estará contenida dentro del bucle while de manera que podrá ser modificada en línea durante el transcurso de la simulación. Figura 3-6. Block Diagram del Control Predictivo Basado en Modelo 35 35 Implementación de controladores predictivos en LabVIEW Para simular la planta de los 4 tanques se llama a través de un bloque Call Library Function Node a la función PLANTA dentro de otra estructura While Loop. Esta estructura se ejecutará cada 100 milisegundos siendo éste el tiempo de integración de la planta. La función está contenida en la misma DLL que el controlador predictivo siendo ésta de la forma 𝑃𝐿𝐴𝑁𝑇𝐴(𝑈, 𝑋, 𝑌) Siendo el argumento 𝑈 la actuación calculada por el controlador y el argumento 𝑌 la salida muestreada de la planta. El estado real de la planta 𝑋 no será en ningún momento accesible por parte del controlador debiendo éste ser inicializado correctamente. El tipo de datos de todos los argumentos serán de nuevo de tipo array data pointer con una resolución double de 8 byte. La actuación calculada por el MPC y la salida de la planta serán de nuevo definidas dentro de una estructura Global Variable que convertirá a dichas variables en globales. Figura 3-7. Block Diagram del modelo no lineal de la planta de los 4 tanques Esta simulación programada en LabVIEW ejecutará por tanto una simulación de la planta de los 4 tanques con un controlador predictivo en tiempo real. Mientras que el modelo no lineal de la planta se ejecuta cada 100 milisegundos, el controlador muestreará la salida de la planta cada 5000 milisegundos y calculará la entrada 𝑈 correspondiente para alcanzar la referencia Ref. Esta simulación en LabVIEW presenta una ventaja ya que si el tiempo de integración de la planta y el tiempo de muestreo del controlador se eligen proporcionales a 100 ms y 5000 ms respectivamente los resultados de la simulación no varían. Es decir, si elegimos como tiempo de integración 20 ms y como tiempo de muestreo 1000 ms obtendremos los mismos datos que para la simulación original solo que la ejecución de la misma tarda 5 veces menos. Durante el transcurso de la simulación, esta aceleración de la misma será configurable a través de una variable global. El tiempo de muestreo y el tiempo de integración de la planta se dividen entre esta variable global de manera que esta variable indica cuantas veces es más rápida la simulación con respecto a la simulación original. En el Front Panel de LabVIEW son configurables antes del inicio de la simulación el estado inicial del sistema, la entrada inicial y siendo la salida inicial una entrada necesaria para el correcto funcionamiento de la DLL pero que no produce ninguna alteración en el resultado de la simulación. El estado estimado incicial del observador, la referencia filtrada inicial, la perturbación estimada inicial y la perturbación filtrada inicial son también configurables para su inicialización. La referencia por su parte, será modificable en cualquier momento una vez que la simulación esté en marcha. Anexo A: Scripts de MATLAB 42 y1 = z1 + (t1 - 1)/(t2)*(z1 - z0); incrementoz = abs(z1 - z0); z0 = z1; t1 = t2; end time(k,3) = toc; if abs(z0 - sol1) < error & abs(zk - sol1) < error k = k + 1; end end %% Resultados gráficos text = {'ISTA','FISTA'}; for k=1:2 figure(k) hold on hf = subplot(1,1,1); bar(time(1:repeticiones,1),'b'); bar(time(1:repeticiones,k+1),'r'); legend('QUADPROG',text{k},'Location','best'); xlabel('Número de problema aleatorio'); ylabel('Tiempo de convergencia (s)','Rotation',90); hf.XLim = [0 repeticiones+1]; hf.FontSize = 25; hf.Box = 'on'; grid on end A2. Comparativa de tiempos de convergencia QUADPROG – FISTA con restricciones de igualdad clc; clear; close all; %% Datos iniciales n = 5; TOL = 1e-5; error = 1e-4; repeticiones = 25; time = zeros(repeticiones,2); k=1; while k <= repeticiones %% Generación del problema D = diag(5*rand(n,1)); A = rand(n-1,n); b = zeros(n-1,1); f = 2*rand(n,1) - ones(n,1); LB = -rand(n,1); UB = 2*rand(n,1); edge = (UB - LB)./2; newz = (UB + LB)./2; %% quadprog(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS) tic z0 = quadprog(D,f,[],[],A,b,LB,UB); time(k,1) = toc; %% ALGORITMO FISTA x0 = zeros(n-1,1); y1 = zeros(n-1,1); t1 = 1; fin = 0; Q = inv(A*inv(D)*A'); 43 tic while fin==0 f2 = (f - A'*y1)'; % Problema 1 z = -f2/D; zaux = min(abs(z' - newz),edge).*sign(z' - newz); zyk = zaux + newz; %%%%%%%%%%%%% if abs(A*zyk - b) < TOL fin = 1; else incrementoy = -Q*(A*zyk - b); x1 = y1 + incrementoy; t2 = 0.5*(1 + sqrt(1 + 4*(t1^2))); y1 = x1 + (t1 - 1)/(t2)*(x1 - x0); x0 = x1; t1 = t2; end end time(k,2) = toc; if abs(z0 - zyk) < error k = k + 1; end end %% Figura QUADPROG - FISTA figure(1) hold on hf = subplot(1,1,1); bar(time(1:repeticiones,1),'b'); bar(time(1:repeticiones,2),'r'); legend('QUADPROG','FISTA','Location','best'); xlabel('Número de problema aleatorio'); ylabel('Tiempo de convergencia (s)','Rotation',90); hf.XLim = [0 repeticiones+1]; hf.FontSize = 25; hf.Box = 'on'; grid on A3. Simulación del modelo no lineal de la planta de los 4 tanques clear; clc; close all; % Puntos de funcionamiento Uf1 = [1.7602; 1.8072]; Uf2 = [1.948; 2]; Uf3 = [1.4802; 1.5197]; Xf1 = puntoFunc(Uf1); Xf2 = puntoFunc(Uf2); Xf3 = puntoFunc(Uf3); %% Parámetros modificables n = 4; m = 2; % Tiempo de muestreo y tiempo de integración Tm = 5; h = 0.1; % Punto de linealización Uf = Uf1; Xf = Xf1; % Referencia alcanzable Ref1 = Xf2; Ref2 = Xf3; Anexo A: Scripts de MATLAB 44 % Estado inicial X = Xf; % Entrada inicial U = Uf; % Observador y perturbación inicial observador = X - Xf; perturb = zeros(m,1); dgorro = zeros(m,1); dr = zeros(m,1); rfiltrada = zeros(m,1); % Parámetros de los filtros rho = 0.8; alpha = 0.8; % Constantes modificables para el MPC N = 20; TOL = 1e-7; repeticionesmax = 2000; Q = 1*eye(n); R = 1e-4*eye(m); Q0 = eye(n+m); R0 = 100*eye(m); %% Parametización del modelo [A,B,C,D] = modelo4tanques(Xf,5); r = Ref1(1:m) - Xf(1:m); r2 = Ref2(1:m) - Xf(1:m); % Cálculo de los límites LBx = 0.2 - Xf; UBx = 1.2 - Xf; LBu = 0 - Uf; UBu = 3 - Uf; % Cálculo de los parámetros del observador A0 = eye(n+m); A0(1:n,1:n) = A; B0 = [C,eye(m)]; Lxd = dlqr(A0',B0',Q0,R0)'; %% Generar parámetros para MPC % Generar límites LB = zeros(N*(m+n),1); UB = zeros(N*(m+n),1); H = zeros(N*(m+n)); k=0; while k <= (N-1)*(m+n) for i=1:n LB(k+i,1) = LBx(i,1); UB(k+i,1) = UBx(i,1); H(k+i,k+i) = Q(i,i); end for j=1:m LB(k+i+j,1) = LBu(j,1); UB(k+i+j,1) = UBu(j,1); H(k+i+j,k+i+j) = R(j,j); end k = k + i + j; end % Generar fronteras edge = (UB-LB)./2; newz = (UB+LB)./2; % Generar matriz E E = kron(eye(N),[A B]); j=0; for i=1:n:n*N-n j = j+1; 45 E(i:i+n-1,((j-1)*(n+m)+(m+n+1)):((j-1)*(n+m)+(n+m+1)+n-1)) = -eye(n); end E = [eye(n) zeros(n,n*N-n+m*N); E]; % Generar vector b b = zeros(n*(N+1),1); % Generar Qfista Qfista = inv(E*inv(H)*E'); %% Parámetros para el SSTO L = [A-eye(n) B;C D]\[zeros(n,m);eye(m)]; c = [L;-L]; d = [UBx;UBu;-LBx;-LBu]; %% Inicialización del sistema NO LINEAL X = simula_NL_4Tanques(X,U,Tm,h)'; Y = C*X; x = X - Xf; y = C*x; u = U - Uf; estado(1,:) = X; entrada(1,:) = U; salida(1,:) = Y; %% Bucle for i=1:1000 %% Observador observador2 = A*observador + B*u + Lxd(1:4,:)*(y - C*observador - dgorro); dgorro = dgorro + Lxd(5:6,:)*(y - C*observador - dgorro); observador = observador2; % Filtro del observador dr = (1 - alpha)*dr + alpha*dgorro; %% Actualización de la referencia if i==500 r = r2; end rk = r - dr; %% Filtro del SSTO rfiltrada = (1 - rho)*rfiltrada + rho*rk; % SSTO indice = c*rfiltrada > d; res = d./(c*rfiltrada); landa = min(res(indice)); if size(landa)==[0 1] landa = 1; end resul = L*rfiltrada*landa; xr = resul(1:n); ur = resul(n+1:n+m); %% MPC % Generar vector f f=[]; for k=1:N f = [f;-Q*xr;-R*ur]; end % Actualización del estado actual b((n*N+1):(n*(N+1))) = xr(1:n); Anexo A: Scripts de MATLAB 46 b(1:n) = observador(1:n); % Algoritmo FISTA x0 = zeros(n*(N+1),1); y1 = zeros(n*(N+1),1); t1 = 1; fin = 0; repeticiones = 0; while (fin==0) f2 = f - E'*y1; % Problema 1 z0 = -f2'/H; zaux = min(abs(z0' - newz),edge).*sign(z0' - newz); z0fista = zaux+newz; %%%%%%%%%%%%% repeticiones = repeticiones + 1; if abs(E*z0fista - b) < TOL | repeticiones > repeticionesmax fin = 1; else incrementoy = -Qfista*(E*z0fista - b); x1 = y1 + incrementoy; t2 = 0.5*(1 + sqrt(1 + 4*(t1^2))); y2 = x1 + (t1 - 1)/(t2)*(x1 - x0); y1 = y2; x0 = x1; t1 = t2; end end u = z0fista(n+1:n+m); %% Sistema NO LINEAL U = u + Uf; X = simula_NL_4Tanques(X,U,Tm,h)'; x = X - Xf; y = C*x; % Perturbación a la salida if rem(i,200)==0 perturb = rand(m,1)/20; end y = y + perturb; Y = y + Xf(1:m); %% Toma de datos estado(i+1,:) = X; obs(i,:) = observador + Xf; entrada(i+1,:) = U; salida(i+1,:) = Y; end %% Representación gráfica % Figura 1 figure(1) for k = 1:4 hf = subplot(2,2,k); plot(1:i,estado(1:i,k),'LineWidth',2); hold on; plot(1:i,obs(1:i,k),'LineWidth',2); legend('Estado','Observador') title(sprintf('Estado %d',k)) xlabel('Intervalo de muestreo') ylabel('X_{0}') hf.FontSize = 22; hf.Box = 'on'; grid on end % Figura 2 figure(2) 47 hf = subplot(1,1,1); plot(1:i,salida(1:i,1),'LineWidth',2); hold on; plot(1:i,salida(1:i,2),'LineWidth',2); legend('Salida 1','Salida 2') plot(1:499,Ref1(1)*ones(1,499),'k--',1:499,Ref1(2)*ones(1,499),'k--'); plot(500:1000,Ref2(1)*ones(1,501),'k--',500:1000,Ref2(2)*ones(1,501),'k--'); xlabel('Intervalo de muestreo') ylabel('Y_{0}') hf.FontSize = 25; hf.Box = 'on'; grid on Anexo A: Scripts de MATLAB 48 49 ANEXO B: DLL B1. Generador de la cabecera %% Genera los datos para la DLL clear; clc; fid = fopen('generar_datos_c.txt','w'); % Puntos de funcionamiento Uf1 = [1.7602; 1.8072]; Uf2 = [1.948; 2]; Uf3 = [1.4802; 1.5197]; Xf1 = puntoFunc(Uf1); Xf2 = puntoFunc(Uf2); Xf3 = puntoFunc(Uf3); %% Parámetros modificables n = 4; m = 2; h = 0.1; % Tiempo de integración en segundos Tm = 5; % Tiempo de muestreo en segundos % Punto de linealización Uf = Uf1; Xf = Xf1; % Parámetros de los filtros rho = 0.7; alpha = 0.8; % Constantes modificables para el MPC N = 20; TOL = 1e-7; repeticionesmax = 2000; % Matrices Q = 5*ones(n,1); R = 3*ones(m,1); Q0 = eye(n+m); R0 = 0.1*eye(m); %% Parametización del modelo [A,B,C,D] = modelo4tanques(Xf,Tm); % Cálculo de los límites LBx = 0.2 - Xf; UBx = 1.2 - Xf; LBu = 0 - Uf; UBu = 3 - Uf; % Generar límites LB = [LBx;LBu]; UB = [UBx;UBu]; % Generar fronteras edge = (UB - LB)./2; newz = (UB + LB)./2; % Cálculo de los parámetros del observador A0 = eye(n+m); A0(1:n,1:n) = A; B0 = [C,eye(m)]; Lxd = dlqr(A0',B0',Q0,R0)'; % Generar matriz H H = zeros(N*(m+n)); k=0; while k <= (N-1)*(m+n) for i=1:n H(k+i,k+i) = Q(i); end for j=1:m H(k+i+j,k+i+j) = R(j); Anexo B: DLL 50 end k = k + i + j; end % Generar matriz E E = kron(eye(N),[A B]); j = 0; for i=1:n:n*N-n j = j+1; E(i:i+n-1,((j-1)*(n+m)+(m+n+1)):((j-1)*(n+m)+(n+m+1)+n-1)) = -eye(n); end E = [eye(n) zeros(n,n*N-n+m*N); E]; % Generar E1 y E2 E1 = E(n+1:2*n,1:2*n+m); ET1 = E(1:2*n,1:n+m)'; ET2 = E(n+1:3*n,n+m+1:2*(n+m))'; % Generar Qfista Qfista = inv(E*inv(H)*E'); %% SSTO L = [A-eye(n) B;C D]\[zeros(n,m);eye(m)]; c = [L;-L]; d = [UBx;UBu;-LBx;-LBu]; %% Generar archivo fprintf(fid,'/* Cabecera fuera de la función */\n'); % Generar dimensiones fprintf(fid,'#define n %d\n',n); fprintf(fid,'#define m %d\n',m); fprintf(fid,'#define N %d\n\n',N); fprintf(fid,'/* Cabecera dentro de la función */\n'); % Generar Uf [ii,~]=size(Uf); fprintf(fid,'double Uf[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',Uf(i)); end fprintf(fid,'%.4f};\n',Uf(ii)); % Generar Xf [ii,~]=size(Xf); fprintf(fid,'double Xf[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',Xf(i)); end fprintf(fid,'%.4f};\n\n',Xf(ii)); % Generar A [ii,jj]=size(A); fprintf(fid,'double A[%d][%d]={{',ii,jj); for i=1:ii for j=1:jj-1 fprintf(fid,'%.4f,',A(i,j)); end if i==ii fprintf(fid,'%.4f}};\n\n',A(ii,jj)); else fprintf(fid,'%.4f},\n\t{',A(i,jj)); end end % Generar B [ii,jj]=size(B); fprintf(fid,'double B[%d][%d]={{',ii,jj); for i=1:ii for j=1:jj-1 fprintf(fid,'%.4f,',B(i,j)); end 51 if i==ii fprintf(fid,'%.4f}};\n\n',B(ii,jj)); else fprintf(fid,'%.4f},\n\t{',B(i,jj)); end end % Generar C [ii,jj]=size(C); fprintf(fid,'double C[%d][%d]={{',ii,jj); for i=1:ii for j=1:jj-1 fprintf(fid,'%.4f,',C(i,j)); end if i==ii fprintf(fid,'%.4f}};\n\n',C(ii,jj)); else fprintf(fid,'%.4f},\n\t{',C(i,jj)); end end % Generar edge [ii,~]=size(edge); fprintf(fid,'double edge[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',edge(i)); end fprintf(fid,'%.4f};\n',edge(ii)); % Generar newz [ii,~]=size(newz); fprintf(fid,'double newz[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',newz(i)); end fprintf(fid,'%.4f};\n\n',newz(ii)); % Generar Q [ii,~]=size(Q); fprintf(fid,'double Q[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',Q(i)); end fprintf(fid,'%.4f};\n',Q(ii)); % Generar R [ii,~]=size(R); fprintf(fid,'double R[%d]={',ii); for i=1:ii-1 fprintf(fid,'%.4f,',R(i)); end fprintf(fid,'%.4f};\n\n',R(ii)); % Generar E1 y E2 [ii,jj]=size(E1); fprintf(fid,'double E[%d][%d]={{',ii,jj); for i=1:ii for j=1:jj-1 fprintf(fid,'%.4f,',E1(i,j)); end if i==ii fprintf(fid,'%.4f}};\n',E1(ii,jj)); else fprintf(fid,'%.4f},\n\t{',E1(i,jj)); end end [ii,jj]=size(ET1); fprintf(fid,'double ET1[%d][%d]={{',ii,jj); for i=1:ii for j=1:jj-1 fprintf(fid,'%.4f,',ET1(i,j)); Anexo B: DLL 58 if (res[i]<TOL) { tol++; } } } /* if abs(E*z0fista-b)<TOL | repeticiones>repeticionesmax */ if ((tol==n*(N+1)) || (repeticiones>repeticionesmax)) { fin=1; } else { /* incrementoy=-Qfista*(E*z0fista-b) */ for(i=0;i<n*(N+1);i++) { /* 84x84-84x1 */ incrementoy[i]=0; for(j=0;j<n*(N+1);j++) { incrementoy[i]+=-Qfista[i][j]*res[j]; } /* x1=y1+incrementoy */ x1[i]=y1[i]+incrementoy[i]; } /* t2=0.5*(1+sqrt(1+4*(t1^2))) */ t2=0.5*(1+sqrt(1+4*t1*t1)); /* y1=x1+(t1-1)/(t2)*(x1-x0) */ for(i=0;i<n*(N+1);i++) { y1[i]=x1[i]+((t1-1)/t2)*(x1[i]-x0[i]); x0[i]=x1[i]; } t1=t2; } } for(i=0;i<m;i++) { U[i]=z0fista[n+i]+Uf[i]; } return 1; } DLLIMPORT double PLANTA (double *U, double *X, double *Y) { int i,j; double h = 0.1; /* Sistema NO LINEAL (LEER Y)*/ X[0] = X[0] + h*((-1.3104e-4*sqrt(2*9.81*X[0])+9.2673e- 5*sqrt(2*9.81*X[2])+(0.3/3600)*U[0])/0.03); X[1] = X[1] + h*((-1.5074e-4*sqrt(2*9.81*X[1])+8.8164e- 5*sqrt(2*9.81*X[3])+(0.4/3600)*U[1])/0.03); X[2] = X[2] + h*((-9.2673e-5*sqrt(2*9.81*X[2])+(0.6/3600)*U[1])/0.03); X[3] = X[3] + h*((-8.8164e-5*sqrt(2*9.81*X[3])+(0.7/3600)*U[0])/0.03); for(j=0;j<4;j++) { if(X[j]<0) { X[j]=0; } } for(i=0;i<m;i++) { Y[i]=X[i]; 59 } return 1; } BOOL APIENTRY DllMain (HINSTANCE hInst /* Library instance handle. */ , DWORD reason /* Reason this function is being called. */ , LPVOID reserved /* Not used. */ ) { switch (reason) { case DLL_PROCESS_ATTACH: break; case DLL_PROCESS_DETACH: break; case DLL_THREAD_ATTACH: break; case DLL_THREAD_DETACH: break; } /* Returns TRUE on success, FALSE on failure */ return TRUE; }