Implementación y validación de controladores predictivos en entornos Siemens

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Implementación y validación de controladores predictivos en entornos Siemens Autor: Francisco Fernández-Llebrez Acedo Tutor: Daniel Limón Marruedo Dpto. Ingeniería Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Implementación y validación de controladores predictivos en entornos Siemens Autor: Francisco Fernández-Llebrez Acedo Tutor: Daniel Limón Marruedo Catedrático Dpto. Ingeniería Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado: Implementación y validación de controladores predictivos en entornos Siemens Autor: Francisco Fernández-Llebrez Acedo Tutor: Daniel Limón Marruedo El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos E ste proyecto ha sido posible gracias a la inconmensurable ayuda brindada por mi tutor, Daniel Limón, y a la desinteresada colaboración del doctorando Víctor Manuel Gracia, quienes me han guiado y respondido a todas mis cuestiones a lo largo de este curso académico. Estoy seguro de que os aguardan grandes éxitos gracias a la pasión y dedicación que ponéis en el mundo de la investigación. A Alicia, a mi familia y a mis amigos, por haberme apoyado y haber confiado en mí durante los momentos más complicados y estresantes de la carrera, que no han sido pocos, y haber celebrado conmigo los logros obtenidos. Gracias de todo corazón, estoy seguro que la vida os brindará grandes experiencias, salud y alegría. Francisco Fernández-Llebrez Acedo Sevilla, 2025 I Resumen E l principal objetivo de este proyecto será la introducción de los conceptos del Control Predictivo Basado en Modelo (MPC) en Controladores Lógicos Programables (PLC), que se presentan como la mejor solución ante la automatización de procesos industriales. En relación a dicho propósito, se hará uso tanto del software de Siemens ® - para la implementación del regulador MPC con TIA PORTAL y la simulación del PLC con PLCSIM Advanced - como de MATLAB ® , donde guardaremos y desde donde envíaremos la evolución de los estados de la planta estudiada, la planta de los cuatro tanques. Todo esto se realizará mediante un servidor OPC UA, que actuará como puente de intercambio de variables en tiempo real. Por último, se hará uso de WinCC, la herramienta de diseño de SCADAs e interfaces HMI de Siemens ® , con un modelo de la planta más realista para poder acercarnos lo máximo posible a la realidad y disponer de un simulador fiel del sistema. III XÍndice Abreviado Apéndice B Códigos TIA PORTAL 59 B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 59 Índice de Tablas 71 Índice de Códigos 73 Bibliografía 75 Glosario 77 Índice Resumen III Abstract V Índice de Figuras VII Índice Abreviado IX Notación XIII 1 Planteamiento y objetivo del proyecto 1 1.1 Principios de la teoría de control automático en variables de estado 3 1.2 Introducción al regulador MPC 4 1.3 La planta de los cuatro tanques 5 2 Control LQR de la planta de los cuatro tanques 9 2.1 El regulador LQR 9 2.1.1 Controlabilidad e implementación de un observador para el regulador LQR 10 2.2 Comunicación MATLAB-PLC vía servidor OPC UA 12 2.3 Implementación del regulador LQR 12 3 Control MPC y simulación de la planta de los cuatro tanques linealizada 19 3.1 El controlador MPC para la planta de los cuatro tanques 19 3.2 El algoritmo ADMM y SPCIES 22 3.2.1 Uso de SPCIES en MATLAB®24 3.3 Control de un PLC mediante un regulador MPC generado con SPCIES 27 3.3.1 Código C para MPC con SPCIES 27 3.3.2 Traducción a SCL 27 3.3.3 Control MPC con PLC simulado 28 3.4 Simulación de la planta ideal 32 4 Control MPC y simulación de la planta no lineal 35 4.1 El modelo no lineal simplificado 35 4.2 Aproximación a la planta real 39 4.3 Implementación de los PID locales en las válvulas 40 4.3.1 Sintonización de los PIDs 41 4.3.2 Envío de las referencias a los PIDs 42 4.3.3 Implementación de los controladores PID 44 4.4 Método alternativo: Sincronización por eventos 45 XI XII Índice 4.5 Simulación de la planta no lineal real 47 5 Conclusiones y líneas futuras 49 Apéndice A Códigos MATLAB®51 A.1 Código control LQR con observador y cancelación de perturbaciones 51 A.2 Simulador de la planta real de los cuatro tanques 55 Apéndice B Códigos TIA PORTAL 59 B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 59 Índice de Tablas 71 Índice de Códigos 73 Bibliografía 75 Glosario 77 Notación RCuerpo de los números reales MPC Model Predictive Control o Control Predictivo Basado en Modelo PLC Programmable Logic Controller o Controlador Lógico Programable LQR Linear Quadratic Regulator o Regulador Cuadrático Lineal SISO Single Input Single Output MIMO Multiple Input Multiple Output ∥v∥Norma del vector v ⟨v,w⟩Producto escalar de los vectores vyw |A|Determinante de la matriz cuadrada A det(A)Determinante de la matriz (cuadrada) A A⊤Transpuesto de A A−1Inversa de la matriz A A†Matriz pseudoinversa de la matriz A c.t.p. En casi todos los puntos c.q.d. Como queríamos demostrar e.o.c. En cualquier otro caso enúmero e ∂y ∂xDerivada parcial de yrespecto a x x◦Notación de grado, xgrados. InMatriz identidad de dimensión n diag(x)Matriz diagonal a partir del vector x diag(A)Vector diagonal de la matriz A :Tal que ∥x∥Norma-2 del vector x xi,i=1,2,...,nElementos i,de1an, del vector x dxDiferencial de x ⩽Menor o igual ⩾Mayor o igual \Backslash ⇔Si y sólo si x=a+3= ↑ a=1 4Igual con explicación XIII XIV Notación a bFracción, a/b ∆Incremento b·10aFormato científico ∈Perteneciente a 1 Planteamiento y objetivo del proyecto "Los grandes resultados requieren grandes ambiciones" — Heráclito D esde los comienzos de la primera revolución industrial, la automatización de las diversas tareas presentadas en las fábricas ha sido una prioridad en la industria. No solo se ha buscado el aumento de la producción y el beneficio, si no también la mejora de las condiciones laborales y la seguridad de los trabajadores, por lo que es importante contar con sistemas de automatización y control robustos y precisos. El tema del que nos ocuparemos durante este proyecto, la búsqueda de soluciones óptimas para el control de automatismos mediante PLCs, se remonta a los finales de la década de 1950 cuando Richard Bellman formula el principio de optimalidad dinámica y a los comienzos de la década de 1960 cuando Rudolf Kálmán desarrolla la teoría del control óptimo lineal cuadrático (LQR). Este regulador se encarga de minimizar una función cuadrática objetivo, conocida como función de coste, cuyos elementos serán los estados y entradas del sistema que estamos estudiando, ponderando los primeros con una matriz Q y las segundas con una matriz R para equilibrar esfuerzo y desempeño del controlador. Aunque el control LQR nos brinda soluciones óptimas para cada tiempo de muestreo, en la industria presenta una gran limitación: no admite restricciones. Por ejemplo, si tuviéramos un tanque lleno de líquido corrosivo que quisiéramos regular hacia un nivel más alto, un LQR no tendría en cuenta restricciones explícitas para llevar el líquido hacia el nivel superior, pudiendo desbordar si en el camino de soluciones óptimas hubiera algunas que superasen la altura del depósito. Debido a esto, uno de los métodos más empleados en la industria actualmente es el MPC o Model Predictive Controller (Controladores Predictivos basados en Modelo). Estos se basan en un modelo del sistema a controlar, de manera que se obtendrán soluciones óptimas de una función de coste similar a la del LQR, pero en un horizonte de tiempo finito y teniendo en cuenta las restricciones deseadas. Así, podremos tener en cuenta limitaciones tanto en las entradas como en los estados.. El LQR y el MPC presentan ciertas ventajas sobre los controladores más empleados, como el PID, especialmente en procesos industriales de mayor complejidad. El LQR es eficaz en la optimización del rendimiento en sistemas lineales y capacidad de ajustar el sistema a condiciones variables aunque el desempeño en tiempo real es limitado por su incapacidad de manejar las restricciones deseadas. El control MPC va un paso más allá para desenvolverse bien en el tiempo real, consideran1 2Capítulo 1. Planteamiento y objetivo del proyecto do como dijimos las restricciones de los estados y las acciones de control en cada tiempo de muestreo. Una de las principales ventajas del MPC es su habilidad para gestionar sistemas multivariables, donde interactúan múltiples entradas y salidas. Este tipo de sistemas es común en la industria moderna, en aplicaciones como el control de tanques, reactores químicos, o procesos de fabricación donde una única variable de control no es suficiente. El MPC puede ajustar simultáneamente múltiples entradas para mantener múltiples salidas dentro de los rangos deseados, optimizando el rendimiento global sin dañar ningún aspecto del sistema. Esto es especialmente crucial en aplicaciones industriales donde los cambios en una variable pueden afectar de forma compleja a otras. Además, el MPC puede trabajar con un horizonte de predicción, lo que significa que no solo optimiza la acción de control en el momento presente, sino que tiene en cuenta las condiciones futuras para tomar decisiones más informadas. Esto es particularmente importante en sistemas donde las dinámicas del proceso son lentas y las decisiones a corto plazo pueden no ser suficientes para mantener el rendimiento a largo plazo. Otra gran ventaja del MPC es su capacidad para incorporar restricciones operacionales. En entornos industriales, las restricciones físicas son una preocupación constante: no se puede permitir que los niveles de un tanque superen un cierto umbral, o que las velocidades de un motor se excedan. El MPC integra estas limitaciones directamente en su función de optimización, asegurando que las soluciones obtenidas sean viables y seguras. Por ejemplo, en el proceso de llenado de un tanque, el MPC puede garantizar que el nivel del líquido nunca supere el 80% de la capacidad del tanque en el camino de soluciones óptimas calculadas, manteniendo la seguridad del proceso sin comprometer el rendimiento. Así, este proyecto no solo busca validar la implementación de MPC en entornos industriales, sino también demostrar que los controladores avanzados como el LQR y el MPC son esenciales para optimizar la automatización de procesos, maximizando el rendimiento, respetando las restricciones operacionales, y ofreciendo soluciones precisas y seguras para los desafíos industriales actuales. Claramente, el hecho de resolver un problema de optimización con restricciones conlleva un costo computacional muy alto, por lo que uno de los problemas a los que se enfrenta la industria es buscar la implementación de controladores MPC cada vez más precisos en dispositivos cada vez más pequeños. En este proyecto observaremos el efecto de los controladores MPC en el entorno de los PLCs de Siemens ® para ver como responderían en la industria, utilizando para ello el programa de configuración de PLCs conocido como TIA PORTAL, y por otra parte el programa de simulación y supervisión de sistemas llamado WinCC. Para ello, emplearemos el software de SPCIES, desarrollado por el GEPOC de la Universidad de Sevilla, por un equipo formado por Pablo Krupa, Víctor Gracia, Daniel Limón y Teodoro Álamo. Su objetivo es el desarrollo de código de controladores MPC para PLCs. Genera código en MATLAB ® y en C, que puede ser adaptado o traducido al lenguaje que necesite el PLC adquirido, que en el caso que nos atañe, será el lenguaje estructurado de Siemens®SCL. Por último, cabe resaltar que emplearemos el entorno de desarrollo de Siemens ® ya que nos ofrece una integración compacta entre todos los programas que vamos a emplear durante el proyecto como TIA PORTAL y PLCSIM, además de una conexión con MATLAB ® a través del soporte OPC UA que explicaremos con detalle más adelante. Esta empresa posee una larga trayectoria como proveedora de soluciones industriales, por lo que sus productos nos garantizan la intercompatibilidad que necesitamos en este proyecto. 1.1 Principios de la teoría de control automático en variables de estado 3 1.1 Principios de la teoría de control automático en variables de estado El control automático es la ciencia que se encarga del modelado de sistemas que regulan de manera automática el funcionamiento de máquinas o procesos. Para ello, partiremos de las ecuaciones diferenciales que caracterizan la evolución del proceso durante el tiempo, conocido como el modelo dinámico de la planta o sistema: ˙x(t) = f(x(t),u(t),t) y(t) = h(x(t),u(t),t) x(t0) = x0 (1.1) Siendo t el tiempo, x(t)∈Rn el vector de estados del sistema, ˙x(t)∈Rn el vector del estado siguiente del sistema, u(t)∈Rmel vector de entradas del sistema, e y(t)∈Rpel vector de salidas. Cabe señalar que en la realidad, los modelos obtenidos con las ecuaciones diferenciales serán no lineales. Sin embargo, una práctica común es simplificar el problema mediante la linealización en torno a un punto de equilibrio. Se podrá representar entonces el sistema en las conocidas ecuaciones de espacio de estados: ˙x(t) = Acx+Bcu y(t) = Ccx+Dcu(1.2) Donde Ac∈Rnxn será la matriz que relaciona el estado siguiente con el estado actual, Bc∈Rnxm la matriz que relaciona la salida con el estado siguiente, Cc∈Rmxn la matriz que relaciona el estado con la salida y Dc∈Rmxm la matriz que relaciona la entrada con la salida del sistema. Para obtener esta ecuación, como mencionamos antes, habrá que definir un punto de funcionamiento o equilibrio sobre el que trabajar. xyuserán las variables absolutas, y xe y ue los puntos de equilibrio: ˜x=x−xe,˜u=u−ue Para la linealización, debemos realizar una expansión en series de Taylor alrededor del punto de equilibrio y deshacernos de los términos no lineales. No nos harán falta términos de segundo orden ya que los de primer orden consiguen una aproximación aceptable: ˙x≈f(xe,ue)+ ∂f ∂x(xe,ue)(x−xe)+ ∂f ∂u(xe,ue)(u−ue) Donde: Ac=∂f ∂x(xe,ue),Bc=∂f ∂u(xe,ue) Y al cumplirse en el equilibrio f(xe,ue) = 0, nos quedará: ˙ ˜x=Ac˜x+Bc˜u De la misma forma para la salida: y≈h(xe,ue)+ ∂h ∂x(xe,ue)(x−xe)+ ∂h ∂u(xe,ue)(u−ue) Donde: Cc=∂h ∂x(xe,ue),Dc=∂h ∂u(xe,ue) 4Capítulo 1. Planteamiento y objetivo del proyecto Por lo que quedará: ˜y=Cc˜x+Dc˜u Emplearemos igualmente la notación de la ecuación (1.2) con x y u como variables incrementales, siendo ˜x y ˜u usadas solo como ejemplo. Observando la ecuación (1.2) , podemos notar que la expresamos en un tiempo continuo. En la realidad, los controladores no pueden muestrear el tiempo de manera continua, por lo que debemos recurrir a métodos de discretización conocidos como el de Euler o el de Tustin. Así, la ecuación en espacio de estados discreta quedará como sigue: x(k+1) = Ax(k)+Bu(k) y(k) = Cx(k)+ Du(k) x(0) = x0 (1.3) Para k ∈Z+ , i.e., k será un número entero no negativo. Ahora los instantes en los que se actualizarán las ecuaciones de estado vendrán dados por t=k·Tm , donde Tm es el tiempo de muestreo seleccionado. Por ejemplo, si Tm = 1 segundo, los estados se actualizarán con una frecuencia de 60 veces por minuto. Para simplificar la notación de la ecuación del espacio de estados, la expresamos de la siguiente manera: x+=Ax +Bu y=Cx +Du x(0) = x0 (1.4) Podemos representar restricciones en nuestros sistemas tanto en las entradas como en las salidas, al poseer los sistemas límites físicos como los niveles máximos de un tanque, la máxima velocidad de un motor o las capacidades límite de los actuadores de nuestra planta. Podemos expresar conjuntos de restricciones como: X={x∈Rn:Axx≤bx} U={u∈Rm:Auu≤bu} donde X y U son los conjuntos de restricciones en los estados y entradas, respectivamente. En el siguiente apartado las particularizaremos para el controlador fundamental de nuestro proyecto, el MPC. Estas restricciones, que incluyen límites físicos como los niveles máximos de un tanque o las capacidades máximas de los actuadores, deben ser consideradas dentro del modelo desde el inicio. De esta forma, el proceso de optimización del MPC no solo busca minimizar una función de coste, sino que también se asegura de que las soluciones obtenidas sean viables dentro de los límites operacionales del sistema, garantizando seguridad y eficiencia en la operación del proceso. 1.2 Introducción al regulador MPC En términos generales, la formulación del problema de un MPC implica: • La predicción de la evolución del sistema utilizando un modelo matemático como el de la ecuación (1.4). • La minimización de una función de coste en un horizonte de predicción finito sujeta a restricciones. J= N−1 ∑ k=0∥x(k)∥2 Q+∥u(k)∥2 R+∥x(N)∥2 P(1.5) 1.3 La planta de los cuatro tanques 5 Siendo Nel horizonte de predicción, x(k) yu(k) los estados y las acciones de control respectivamente a lo largo de Nyx(N) el estado terminal. • La aplicación solo de la primera acción de control, repitiendo el proceso en el siguiente instante de muestreo. Las matrices de la función de coste (1.5) , Q y R, positivas semidefinidas, servirán para la ponderación de las desviaciones con respecto a la referencia y la penalización del esfuerzo de control respectivamente, y la matriz P, será la matriz terminal que penaliza el estado final del horizonte. Esta penalización terminal con P mejora la estabilidad del sistema y es, en muchos casos, elegida como la solución de una ecuación de Riccati, al igual que en el controlador LQR, el cual explicaremos con detalle en el siguiente capítulo. Los controladores MPC calcularán las predicciones partiendo de la expresión de la ecuación (1.4) , que indica el valor del estado siguiente con respecto a estados y entradas anteriores. Siguiendo esta recursión: x(1) = Ax(0)+Bu(0) x(2) = A2x(0)+ABu(0)+Bu(1) x(3) = A3x(0)+A2Bu(0)+ABu(1)+Bu(2) . . . x(N) = ANx(0)+AN−1Bu(0)+AN−2Bu(1)+···+Bu(N−1)(1.6) Que se podrá expresar matricialmente como sigue:        x(1) x(2) x(3) . . . x(N)        | {z } xF =       B0 0 ... 0 AB B 0... 0 A2B AB B ... 0 . . .. . .. . ..... . . AN−1B AN−2B AN−3B... B        | {z } Gu ∗       u(0) u(1) u(2) . . . u(N−1)        | {z } uF +       A A2 A3 . . . AN        |{z} Gx x(0)(1.7) xF=GuuF+Gxx(0)(1.8) Las restricciones a las que estará sujeta la función de coste (1.5) se pueden denotar como sigue: x(k)∈X≜{x∈Rn:Axx≤bx} u(k)∈U≜{u∈Rm:Auu≤bu}(1.9) En la que nos detendremos luego con más detalle y particularizaremos para nuestra planta, que introduciremos a continuación. 1.3 La planta de los cuatro tanques Esta planta, la planta de los cuatro tanques, será la empleada para estudiar los efectos de los distintos controladores que tratará este proyecto ya que, al ser una planta multivariable con dos entradas y cuatro salidas, posee un notable interés como problema de control de complejidad media. Está compuesta por dos tanques inferiores (1 y 2), y otros dos superiores (3 y 4) que descargan sobre los inferiores. Los tanques se llenarán a través de los caudales qa y qb , que como se puede apreciar, se dividirán en dos ramas. Así, el volumen del tanque 1 crecerá como qaγa , el tanque 3 como qa(1−γa) , 12 Capítulo 2. Control LQR de la planta de los cuatro tanques la convergencia del error a 0 será por ubicación de polos, escogiendo polos estables (dentro del círculo unidad del plano complejo) y más rápidos que los del sistema. También se podría calcular la ganancia óptima de manera similar a como lo hicimos con el regulador LQR, resolviendo la ecuación de Ricatti. Se adjuntará en el apéndice un ejemplo de control LQR más observador por ubicación de polos con cancelación de perturbaciones para la planta de los cuatro tanques. 2.2 Comunicación MATLAB-PLC vía servidor OPC UA Para llevar a cabo la comunicación entre un PLC como actuador y MATLAB como simulador, emplearemos un servidor OPC UA (Unified Architecture). Estos han sido muy importantes en la automatización industrial reciente al permitir la comunicación entre distintos tipos de dispositivos, sin importar la plataforma o el fabricante. El entorno de programación Siemens ® que emplearemos para esta prueba será TIA PORTAL (Totally Integrated Automation Portal), en el que podremos configurar tanto el hardware como software de nuestro PLC. Cabe destacar que en primer lugar se pretendía emplear TIA PORTAL V13 con un PLC S7-1200, aunque pudimos observar que tanto esta versión de TIA PORTAL como esta gama de PLCs no incluían de ninguna manera las opciones OPC UA que nos interesaban. A modo de resumen de las compatibilidades se presentan estas tablas: Tabla 2.1 Compatibilidad de PLC Siemens con OPC UA. Modelo de PLC Soporte OPC UA S7-1500 Sí (puede actuar como servidor OPC UA) S7-1200 No (no puede ser servidor OPC UA) Otros modelos No (requiere gateway o hardware adicional) Tabla 2.2 Soporte de OPC UA según versión de TIA Portal. Versión de TIA Portal Soporte OPC UA V13 No (no incluye opciones OPC UA) V14–V16 Parcial (requiere licencias adicionales) V17 y posteriores Sí (soporte nativo y completo de OPC UA) Por ello, la versión de TIA PORTAL seleccionada finalmente será la V17 y la CPU usada para las pruebas será la 1516-3 PN/DP, pertenenciente a la gama de los PLC S7-1500 de Siemens ® , cuya versión firmware (V 2.9) nos permite la opción de habilitar al mismo PLC como servidor OPC UA. Así, empleando como simulador de PLC el programa PLCSIM Advanced, también de Siemens ® , podremos establecer una comunicación entre un PLC simulado que hará de servidor OPC UA y MATLAB ® , que adoptará el papel de cliente OPC UA gracias al Toolbox de comunicaciones industriales Industrial Communication. 2.3 Implementación del regulador LQR Iniciados ya en la teoría del controlador LQR y en el método elegido para la comunicación entre PLC y MATLAB ® , quedará realizar los códigos en TIA PORTAL y MATLAB ® . El esquema que vamos a seguir será del tipo Hardware in the Loop, en el que validaremos un controlador implementado sobre un hardware, en este caso el PLC, sobre una planta simulada. A continuación, se muestra un esquema simplificado: 2.3 Implementación del regulador LQR 13 Figura 2.2 Esquema de conexión entre PLC (servidor OPC UA) y MATLAB. En primer lugar, en TIA PORTAL habrá que realizar la declaración de variables y publicarlas en el servidor OPC UA para que el cliente MATLAB pueda acceder a ellas. Como queremos seguir el control por realimentación de estados de la ecuación (2.3) , las variables que emplearemos serán un vector de estados actuales, otro para el estado y entrada de equilibrio y otro para las entradas o acciones de control calculadas. Al ser el PLC el actuador, este se encargará de calcular la acción de control y, al ser MATLAB ® el simulador, deberá ir notificando las actualizaciones en el estado actual. La sincronización entre MATLAB ® y el PLC no se realiza de manera continua en el tiempo, sino que es por eventos. Esto significa que no hay una sincronización estricta en tiempo entre el modelo y el controlador. En lugar de un muestreo regular, el PLC realiza sus cálculos de control y envía las actualizaciones de los estados al cliente OPC UA cuando ha terminado de calcular una nueva acción de control. Este método de sincronización por eventos tiene la ventaja de evitar una simulación en tiempo real, permitiendo realizar nuestros experimentos en un tiempo acortado y garantizando que las acciones solo se envíen cuando se han calculado nuevas soluciones de control. Para gestionar esta sincronización, se utiliza una variable auxiliar que desempeña la función de notificar el fin de cálculo. Esta variable le indica a MATLAB ® cuándo el PLC ha terminado de calcular la acción de control y está listo para recibir el siguiente conjunto de actualizaciones de estado. Sin esta variable de notificación, dado que MATLAB ® tiene una capacidad de procesamiento mucho mayor que el PLC, podrían enviarse múltiples estados durante un solo ciclo de control, lo que generaría una desincronización y posibles pérdidas de datos durante la simulación. Al sincronizar ambos sistemas con esta señal de notificación, se asegura que MATLAB ® solo envíe los estados cuando el PLC ha completado el cálculo, evitando así cualquier sobresaturación de datos. El tiempo de muestreo del sistema se fija en MATLAB ® y se configura como un parámetro dentro del modelo. Aunque el tiempo de muestreo será independiente del ciclo de control del PLC, habrá que asegurar que ambos sistemas trabajen de manera coordinada para que el modelo pueda procesar de manera ordenada las acciones de control que le llegan. Es importante además conocer que las tasas de escritura y lectura del servidor OPC UA serán configurables, aunque normalmente los tiempos se configurarán en torno a los 0.5 segundos. El código en TIA Portal será bastante simple, al poder además calcular la ganancia K offline con la función dlqr que ya se mencionó, permitiendo al PLC realizar los cálculos de control de forma eficiente sin necesidad de calcular la ganancia en tiempo real. // Calcular u = -K * (x - x_eq) + u_eq IF NOT #accioncalculada THEN 14 Capítulo 2. Control LQR de la planta de los cuatro tanques #accion[0] := #accioneq[0] - #K00 * (#estado[0] - #estadoeq[0]) - #K01 * (# estado[1] - #estadoeq[1]) - #K02 * (#estado[2] - #estadoeq[2]) - #K03 * (#estado[3] - #estadoeq[3]); #accion[1] := #accioneq[1] - #K10 * (#estado[0] - #estadoeq[0]) - #K11 * (# estado[1] - #estadoeq[1]) - #K12 * (#estado[2] - #estadoeq[2]) - #K13 * (#estado[3] - #estadoeq[3]); END_IF; #accioncalculada := TRUE; Código 2.1 Control LQR parte PLC en TIA PORTAL. 2.3 Implementación del regulador LQR 15 Y el bloque principal, con sus entradas y salidas será: Figura 2.3 Bloque principal control LQR. Por otra parte, en MATLAB ® quedará como sigue, suponiendo la planta ya inicializada (Código 1.1). Veremos el código en dos partes: 1% Configuración de conexión OPC UA 2opcUaClient = opcua(’opc.tcp://192.168.0.2:4840’); % Reemplaza con la IP del PLC 3connect(opcUaClient); 4 5% Variables a las que podemos acceder 6ServerNodes = browseNamespace(opcUaClient); 7disp(ServerNodes) 8 9 10 % Definir nodos OPC UA 11 nodoEstado = findNodeByName(ServerNodes, ’estado’); 12 nodosEstado = nodoEstado.Children; % Extraer los subnodos 13 14 % Obtener nodos hijos de accioneq 15 nodoAccion = findNodeByName(ServerNodes, ’accion’); 16 nodosAccion = nodoAccion.Children; % Extraer los subnodos 17 18 nodoEstadoEq = findNodeByName(ServerNodes, ’estadoeq’); 19 nodosEstadoEq = nodoEstadoEq.Children; % Extraer los subnodos 20 21 % Obtener nodos hijos de accioneq 22 nodoAccionEq = findNodeByName(ServerNodes, ’accioneq’); 23 nodosAccionEq = nodoAccionEq.Children; % Extraer los subnodos 24 25 %Nodo flag accioncalculada 26 nodoAccionCalculada = findNodeByName(ServerNodes, ’accioncalculada’); Código 2.2 Acceso a las variables del servidor OPC UA. En esta primera parte, accedemos con los comandos que se observan a las variables publicadas en el servidor OPC UA, que es el mismo PLC. Tras esto, ejecutaremos el control LQR para verificar la validez de nuestro sistema: 1%Control LQR 16 Capítulo 2. Control LQR de la planta de los cuatro tanques 2 3%matriz de coste del estado Q y matriz coste entrada R 4[K]=dlqr(Ad,Bd,Q,R); 5 6%xeq=[Ad-eye(4),B; Cd, Dd]\[zeros(4, 1);ref]; 7 8% Bucle de control LQR 9for k = 1:Tsim 10 if k >= (Tsim/3) && k < (2*Tsim/3) 11 ref = [1.1; 0.8]; 12 nueq = [Ad-eye(4),Bd; Cd, Dd]\[zeros(4, 1);ref]; 13 xeq = nueq(1:4); 14 ueq = nueq(5:6); 15 for i = 1:4 16 writeValue(opcUaClient, nodosEstadoEq(i), double(xeq(i))); 17 end 18 for i = 1:2 19 writeValue(opcUaClient, nodosAccionEq(i), double(ueq(i))); 20 end 21 elseif k >= (2*Tsim/3) 22 ref = [0.8; 1]; 23 nueq = [Ad-eye(4),Bd; Cd, Dd]\[zeros(4, 1);ref];%calculamos alturas equilibrio y caudales equilibrio para esta referencia 24 xeq = [nueq(1);nueq(2);nueq(3);nueq(4)]; 25 ueq = [nueq(5);nueq(6)]; 26 for i = 1:4 27 writeValue(opcUaClient, nodosEstadoEq(i), double(xeq(i))); 28 end 29 for i = 1:2 30 writeValue(opcUaClient, nodosAccionEq(i), double(ueq(i))); 31 end 32 else 33 ref = [0.5; 0.5]; % Punto de operación 34 end 35 36 37 % Fórmula LQR, a partir de Tsim/3 entra en funcionamiento el control 38 if k >= (Tsim/3) 39 40 41 for i=1:4 42 writeValue(opcUaClient, nodosEstado(i), double(xlqr(i))); 43 end 44 pause(0.1); % Pequeña pausa para no sobrecargar la comunicación 45 46 accioncalculada = readValue(opcUaClient, nodoAccionCalculada); 47 48 if (accioncalculada==1) 49 for i = 1:2 50 uk(i) = readValue(opcUaClient, nodosAccion(i)); 51 end 52 end 53 accioncalculada = 0; 54 writeValue(opcUaClient, nodoAccionCalculada, accioncalculada); 55 56 xlqr = Ad * xlqr + Bd * uk ; 2.3 Implementación del regulador LQR 17 57 58 end 59 60 61 Mxlqr(:, k + 1) = xlqr; %Almacenamiento de históricos 62 Qx(:, k+1) = uk; 63 end Código 2.3 Control LQR parte MATLAB. Podemos observar como en primer lugar nos centramos en la obtención del punto de equilibrio a seguir cada vez que cambia la referencia, y tras esto, aplicaremos la acción de control calculada por el PLC que leeremos en cada tiempo de muestreo para poder obtener el estado siguiente. Ejecutando el script de MATLAB ® y con la simulación activa del PLC, para una relación R/Q aproximadamente igual a 1, los resultados obtenidos quedan como se muestra en estas gráficas: Figura 2.4 Evolución del sistema con regulador LQR. Podemos además comparar los resultados del control MATLAB ® - PLC (izquierda), con tanto la simulación como el control realizados únicamente en MATLAB (derecha): Figura 2.5 Comparación control PLC y control MATLAB. Se puede apreciar como, aún teniendo MATLAB una capacidad de procesamiento mayor, el PLC consigue unos resultados prácticamente iguales para esta planta de los cuatro tanques. El código implementado en el PLC es bastante simple y al ser la ganancia LQR calculada de manera offline, 18 Capítulo 2. Control LQR de la planta de los cuatro tanques podíamos esperar este resultado sin anomalías. Otra observación a realizar, será el incumplimiento de las restricciones para esta referencia ya que podemos notar como los tanques 3 y 4 superarán los niveles permitidos de 1.2 metros de altura al no tener limitaciones explícitamente definidas en nuestro regulador. Trataremos, en el siguiente capítulo, este desajuste con el regulador MPC. 3 Control MPC y simulación de la planta de los cuatro tanques linealizada La medida de la inteligencia es la capacidad de cambiar — Albert Einstein E n este capítulo se tratará la teoría del control MPC aplicada a la planta de los cuatro tanques ideal y linealizada. De la misma forma que en el capítulo anterior, el objetivo será implementar el regulador MPC en un PLC simulado a través de TIA PORTAL, empleando MATLAB ® como simulador de la planta. 3.1 El controlador MPC para la planta de los cuatro tanques En primer lugar, presentaremos el problema de optimización a resolver que vendrá dado por: m´ ın u N−1 ∑ i=0∥x(i)−xr∥2 Q+∥u(i)−ur∥2 R+∥x(N)−xr∥2 P(3.1) s.t. x(i+1) = Ax(i)+Bu(i),i=0,...,N−1(3.2) x(0) = x0(3.3) x(i)∈X,i=0,...,N−1(3.4) u(i)∈U,i=0,...,N−1(3.5) x(N)∈ΩK(3.6) Donde: •x(i)yu(i)son los vectores de estados y entradas en el instante i. •xr(i)yur(i)son las referencias deseadas para los estados y entradas en el instante i. 19 20 Capítulo 3. Control MPC y simulación de la planta de los cuatro tanques linealizada •QyRson las matrices de ponderación descritas en la introducción. •Ppenaliza la diferencia en el estado final x(N)con respecto a la referencia final xr(N). •XyUson los conjuntos de restricciones en los estados y entradas. •ΩKes el conjunto de restricciones sobre el estado final. Este tipo de problemas son denominados QP o Quadratic Programming, y son aquellos en los que la función objetivo es cuadrática, i.e., tiene términos de segundo orden, y las restricciones son lineales. En su forma estándar, un problema de QP se puede expresar como: m´ ın x 1 2xTQx+cTx sujeto a: Ax≤b,x∈X Donde: •xes el vector de variables de decisión. •Q es una matriz simétrica positiva definida que describe la parte cuadrática de la función objetivo. •ces un vector que representa el término lineal de la función objetivo. •Aybrepresentan las restricciones lineales sobre las variables de decisión. Es por esto que la formulación QP es ampliamente utilizada en el campo de los MPC, ya que con esta podremos resolver problemas en tiempo real atados a unas restricciones determinadas. Además, tendremos la capacidad de penalizar como deseemos los estados y entradas de nuestras funciones objetivo mediante las matrices de ponderación explicadas. Revisando la ecuación (1.7) y las matrices A y B de nuestra planta, tendremos la expresión matricial del sistema de los cuatro tanques en un horizonte de predicción N, que define la evolución de los N estados siguientes en función de las entradas y el estado inicial actual:        x(1) x(2) x(3) . . . x(N)        | {z } xF =       B0 0 ... 0 AB B 0... 0 A2B AB B ... 0 . . .. . .. . ..... . . AN−1B AN−2B AN−3B... B        | {z } Gu ∗       u(0) u(1) u(2) . . . u(N−1)        | {z } uF +       A A2 A3 . . . AN        |{z} Gx x(0)(1.7) Con AyBque fueron ya obtenidas a partir de la discretización de AcyBcsiendo estas: Ac=                −a1g Aq2gheq 1 0a3g Aq2gheq 3 0 0−a2g Aq2gheq 2 0a4g Aq2gheq 4 0 0 −a3g Aq2gheq 3 0 0 0 0 −a4g Aq2gheq 4                3.1 El controlador MPC para la planta de los cuatro tanques 21 Bc=          γa 3600A0 0γb 3600A 01−γb 3600A 1−γa 3600A0          Y sus restricciones, descritas en la ecuación (1.11): 0.2m≤hi≤1.2mi=1,2,3,4 0m3 h≤qz≤3m3 hz=a,b(1.11) Como se mencionó antes, las restricciones para los estados y las salidas del horizonte de predicción se pueden describir así: x(k)∈X≜{x∈Rn:Axx≤bx} u(k)∈U≜{u∈Rm:Auu≤bu}(3.7) Donde el vector bx , vector de restricción de estado, estará explícitamente compuesto por las restricciones que deseamos aplicar a nuestro problema, y la matriz Ax se compondrá de dos matrices identidad para expresar la inecuación en forma matricial simplemente. Por ejemplo, para un sistema de cuatro estados: I4 −I4x(k)≤xm´ ax −xm´ ın(3.8)             1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 −1 0 0 0 0−1 0 0 0 0 −1 0 000−1             | {z } Ax     x1 x2 x3 x4    ≤             xm´ ax1 xm´ ax2 xm´ ax3 xm´ ax4 −xm´ ın1 −xm´ ın2 −xm´ ın3 −xm´ ın4             | {z } bx (3.9) Donde nos referimos a x1 , x2 , x3 , x4 como las componentes del vector de estado. Generalizando para el caso del horizonte de predicción igual a N:        Ax0 0 ... 0 0Ax0... 0 0 0 Ax... 0 . . .. . .. . ..... . . 0 0 0 ... Ax        | {z } AxF        x(1) x(2) x(3) . . . x(N)        | {z } xF ≤       bx bx bx . . . bx        |{z} bxF (3.10) Donde ahora nos referimos a x(1) , x(2) , etc., como los vectores de estado en sí, en los periodos 1, 2, 3..., Ndespués del instante actual muestreado. De esta manera, nos aseguramos de imponer los límites deseados para todo el horizonte de predicción, y no solo sobre la acción y estado actuales calculados. 28 Capítulo 3. Control MPC y simulación de la planta de los cuatro tanques linealizada 1// Update the reference 2for(unsigned int j = 0; j < nn_; j++){ 3q[j] = Q[j]*xr[j]; 4qT[j] = 0.0; 5for(unsigned int i = 0; i < nn_; i++){ 6qT[j] = qT[j] + T[j][i]*xr[i]; 7} 8} 9for(unsigned int j = 0; j < mm_; j++){ 10 q[j+nn_] = R[j]*ur[j]; 11 } Código 3.5 Código original generado en C. // Update the reference FOR #jj := 0 TO #n-1DO #q[#jj] := #QQ[#jj] * #xr[#jj]; #qT[#jj] := 0.0; FOR #ii := 0 TO #n-1DO #qT[#jj] := #qT[#jj] + #TT[#jj, #ii] * #xr[#ii]; END_FOR; END_FOR; FOR #jj := 0 TO #m-1DO #q[#jj + #n] := #RR[#jj] * #ur[#jj]; END_FOR; Código 3.6 Código traducido a SCL. Se pueden apreciar a simple vista cambios en la estructura de los bucles for o en las variables, que en lenguaje SCL deberán ir siempre precedidas del símbolo #. 3.3.3 Control MPC con PLC simulado Una vez explicado el proceso para obtener el código en un lenguaje que pueda procesar nuestro PLC, será el momento de regular la altura de los tanques con este. Debido a la extensión del código traducido a SCL, este se adjuntará en el apéndice y se podrá consultar a la vez que se observa el código MATLAB ® . Podremos apreciar las entradas y salidas del bloque de función MPC en esta imagen: Para el código MATLAB ® , seguiremos el mismo procedimiento que con el regulador LQR. En primer lugar, tendremos que acceder a los nodos publicados en el servidor OPC UA del PLC. En esta ocasión, las variables new_data_available y new_state_received, serán las encargadas de coordinar las acciones calculadas por el PLC con los nuevos estados que va calculando MATLAB®. 1clear 2close all 3.3 Control de un PLC mediante un regulador MPC generado con SPCIES 29 Figura 3.3 Bloque de función MPC de TIA PORTAL. 3%Problema de los 4 tanques para qa0=1.948 y qb0=2 4Tm = 5;%tiempo muestreo 5Tsim =5000;%tiempo total de la simulacion 6Nsim = Tsim/Tm; 7 8% Configuración de conexión OPC UA 9opcUaClient = opcua(’opc.tcp://192.168.0.2:4840’); % Reemplaza con la IP del PLC 10 connect(opcUaClient); 11 12 % Variables a las que podemos acceder 13 ServerNodes = browseNamespace(opcUaClient); 14 disp(ServerNodes) 15 16 17 % Definir nodos OPC UA 18 nodoEstado = findNodeByName(ServerNodes, ’estado’); 19 nodosEstado = nodoEstado.Children; % Extraer los subnodos 20 21 % Obtener nodos hijos de accioneq 22 nodoAccion = findNodeByName(ServerNodes, ’accion’); 23 nodosAccion = nodoAccion.Children; % Extraer los subnodos 24 25 nodoEstadoEq = findNodeByName(ServerNodes, ’estadoeq’); 26 nodosEstadoEq = nodoEstadoEq.Children; % Extraer los subnodos 27 28 % Obtener nodos hijos de accioneq 29 nodoAccionEq = findNodeByName(ServerNodes, ’accioneq’); 30 nodosAccionEq = nodoAccionEq.Children; % Extraer los subnodos 31 32 %Nodo flag accioncalculada 33 nodoAccionCalculada = findNodeByName(ServerNodes, ’new_data_available’); 34 35 %Nodo flag estadocalculado 36 nodoEstadoCalculado = findNodeByName(ServerNodes, ’new_state_received’); Código 3.7 Acceso a los nodos publicados de TIA PORTAL. Una vez creadas las matrices A, B, C, D discretas y variables necesarias, como ya se hizo en el capítulo dos, podremos proceder a realizar el bucle de control. Igual que para el control LQR, cambiaremos la referencia en dos instantes concretos. Como podemos notar, cuando escribimos en la variable estadocalculado (asociada a new_state_received) un 1, dará comienzo el funcionamiento del PLC: se restablecerá el valor de las variables fundamentales y el número de iteraciones a 0, y 30 Capítulo 3. Control MPC y simulación de la planta de los cuatro tanques linealizada esperará a que MATLAB ® le envíe un nuevo estado, es decir, que new_state_received = 2. Con el estado recibido, el PLC podrá empezar a calcular la solución subóptima que cumpla con los márgenes de tolerancia como se aprecia en el código SCL completo de TIA PORTAL del apéndice. Una vez tenemos la solución, el PLC cambiará new_data_available de 2 a 1, lo que expulsará a MATLAB ® de la espera infinita del bucle while-end para que pueda calcular un nuevo estado y enviárselo al PLC, cerrando el ciclo. 1%Control MPC 2 3%matriz de coste del estado Q y matriz coste entrada R 4[K]=dlqr(Ad,Bd,Q,R); 5 6estadocalculado = 0; 7writeValue(opcUaClient, nodoEstadoCalculado, estadocalculado); 8 9% Bucle de control MPC 10 for k = 1:Nsim 11 if k >= (Nsim/3) && k < (2*Nsim/3) 12 ref = [0.5; 0.75]; 13 nueq = [Ad-eye(4),Bd; Cd, Dd]\[zeros(4, 1);ref]; 14 xeq = nueq(1:4); 15 ueq = nueq(5:6); 16 for i = 1:4 17 writeValue(opcUaClient, nodosEstadoEq(i), double(xeq(i))); 18 end 19 for i = 1:2 20 writeValue(opcUaClient, nodosAccionEq(i), double(ueq(i))); 21 end 22 23 elseif k >= (2*Nsim/3) 24 ref = [1.1; 0.7]; 25 nueq = [Ad-eye(4),Bd; Cd, Dd]\[zeros(4, 1);ref];%calculamos alturas equilibrio y caudales equilibrio para esta referencia 26 xeq = [nueq(1);nueq(2);nueq(3);nueq(4)]; 27 ueq = [nueq(5);nueq(6)]; 28 for i = 1:4 29 writeValue(opcUaClient, nodosEstadoEq(i), double(xeq(i))); 30 end 31 for i = 1:2 32 writeValue(opcUaClient, nodosAccionEq(i), double(ueq(i))); 33 end 34 else 35 ref = [0.0; 0.0]; 36 end 37 38 39 % A partir de Tsim/3 entra en funcionamiento el control 40 if k >= (Nsim/3) 41 42 estadocalculado = 1; 43 writeValue(opcUaClient, nodoEstadoCalculado, estadocalculado); %asi el PLC comienza a calcular 44 45 for i=1:4 46 writeValue(opcUaClient, nodosEstado(i), double(xlqr(i))); 47 end 3.3 Control de un PLC mediante un regulador MPC generado con SPCIES 31 48 pause(0.1); % Pequeña pausa para no sobrecargar la comunicación 49 50 estadocalculado = 2; %para entrar en la parte iterativa del algoritmo 51 writeValue(opcUaClient, nodoEstadoCalculado, estadocalculado); %asi el PLC comienza a calcular 52 53 while(readValue(opcUaClient, nodoAccionCalculada)==2)%asi mientras newdataavailable=2, esperamos 54 end 55 56 if (readValue(opcUaClient, nodoAccionCalculada)==1) 57 for i = 1:2 58 uk(i) = readValue(opcUaClient, nodosAccion(i)); 59 end 60 end 61 % accioncalculada = 0; 62 % writeValue(opcUaClient, nodoAccionCalculada, accioncalculada); 63 64 xlqr = Ad * xlqr + Bd * uk ; 65 66 end 67 68 Mxlqr(:, k + 1) = xlqr; %Almacenamiento historicos 69 Qx(:, k+1) = uk; 70 end Código 3.8 Bucle de control MPC. 32 Capítulo 3. Control MPC y simulación de la planta de los cuatro tanques linealizada Representando gráficamente los resultados para las referencias señaladas en este experimento, para N=10,ρ=5y una relación R Qsimilar a 1, obtenemos: Figura 3.4 Resultados MPC MATLAB-TIA PORTAL. Resaltando de nuevo en las gráficas como tanto la acción de control como las alturas que podrán alcanzar los tanques estarán restringidos, priorizando las limitaciones antes que el seguimiento de la referencia. Cabe mencionar como en una primera prueba se experimentó con un valor de ρ=200 , lo cual inducía una alta rapidez en la convergencia que desestabilizaba al algoritmo, llevándolo a valores inalcanzables para el PLC y provocando la obtención de Nan (Not a Number) en algunas variables. Esto refuerza la importancia de escoger valores adecuados de los parámetros iniciales o invariables en el tiempo si no queremos que las acciones de control calculadas por el PLC puedan volverse totalmente inestables. Para el ejemplo funcional con ρ=5 , hemos podido apreciar un tiempo de computación más alto debido a la capacidad inferior de procesamiento que posee el PLC y a las pausas realizadas para no saturar la comunicación PLC-MATLAB ® . A pesar de esto, los resultados son bastante satisfactorios y trataremos de realizar una simulación gráfica añadiendo al esquema WinCC, una plataforma de Siemens®que permite la visualización, el control y la supervisión de procesos industriales. 3.4 Simulación de la planta ideal En este apartado final se tratará de realizar una simulación visual de la evolución de las alturas de los tanques con WinCC (Windows Control Center), que es una plataforma de visualización industrial (HMI/SCADA) desarrollada por Siemens ® . Permite el diseño, configuración y supervisión de interfaces hombre-máquina (HMI) para sistemas de automatización. Emplearemos la versión de WinCC unificada con TIA PORTAL, aunque cabe mencionar que existe una versión de WinCC por separado, denominada WinCC Explorer. En primer lugar, se realizará una simulación de la planta teórica, con las válvulas de tres vías y las bombas de caudal qa yqb, descrita en los primeros apartados del proyecto. Como tenemos el diseño de la planta, bastará con seleccionar en WinCC los elementos pertinentes de la librería gráfica de Siemens ® en este caso, ya que es una planta simple y no nos hace falta ningún icono externo. La planta queda como sigue: 3.4 Simulación de la planta ideal 33 Figura 3.5 Planta teórica de los cuatro tanques. Ahora, quedará establecer la conexión de las variables, especialmente las de estado, de nuestro proyecto en el HMI para poder observar la variación de las alturas. Esta conexión se realizará distinguiendo entre variables HMI y variables PLC, y se asignará a las primeras las direcciones deseadas de las segundas mediante direccionamiento simbólico, i.e., las variables HMI accederán por nombre a las del PLC para que en caso de que haya cambio de direcciones, las variables puedan continuar conectadas entre sí. Figura 3.6 Planta teórica de los cuatro tanques simulada. En esta imagen se puede apreciar la planta en simulación para una prueba en la que el tanque uno alcanzará los 0.5 metros y el tanque dos, los 0.75 metros. 4 Control MPC y simulación de la planta no lineal Una vez que aceptamos nuestros límites, los superamos — Albert Einstein L os objetivos de este capítulo final apuntarán a un acercamiento a la realidad, tomando ahora el modelo no lineal de la planta de los cuatro tanques. En primer lugar, trataremos de simular el modelo no lineal ideal, más simplificado, sin la dinámica de las válvulas, para posteriormente tratar de representar lo más fielmente posible la planta de los cuatro tanques ubicada en el laboratorio, que será explicada con detalle en este último capítulo. 4.1 El modelo no lineal simplificado Este modelo de la planta de los cuatro tanques proviene directamente de las ecuaciones que empleamos para la linealización del sistema, que como ya vimos se regían según: Adh1 dt =−a1p2gh1+a3p2gh3+γa qa 3600, Adh2 dt =−a2p2gh2+a4p2gh4+γb qb 3600, Adh3 dt =−a3p2gh3+(1−γb)qb 3600, Adh4 dt =−a4p2gh4+(1−γa)qa 3600. (4.1) Con estas ecuaciones, que no incluyen ni la dinámica de las válvulas, ni las pérdidas de carga en las tuberías, ni la formación de vórtices en los depósitos, podremos implementar en MATLAB ® una función cuyas entradas sean los caudales y alturas actuales, y nos devuelva las alturas para el siguiente instante de muestreo. 1 2function [dz,qv,h]=Modelo4TIdeal(z,q) 3%% Modelo ideal de la planta de los 4 tanques 4 5% Unidades: 6% Alturas: metros 7% Caudales: m3/h 35 36 Capítulo 4. Control MPC y simulación de la planta no lineal 8% Presión: bar 9gama = 0.3; 10 gamb = 0.4; 11 12 %% Cálculo de la Dinámica de la planta 13 14 qv(1) = gama*q(1); 15 qv(2) = gamb*q(2); 16 qv(3) = (1-gamb)*q(2); 17 qv(4) = (1-gama)*q(1); 18 19 %Tanques 20 21 At=0.03; 22 Kt=sqrt(2*9.81)*[1.3104,1.5074,0.92673,0.88164]*1e-4; 23 h=z; 24 h=min(max(h,0.05),1.3); 25 dh=h; 26 27 dh(1)=1/At*(-Kt(1)*sqrt(h(1)) + Kt(3)*sqrt(h(3)) + qv(1)/3600); 28 dh(2)=1/At*(-Kt(2)*sqrt(h(2)) + Kt(4)*sqrt(h(4)) + qv(2)/3600); 29 dh(3)=1/At*(-Kt(3)*sqrt(h(3)) + qv(3)/3600); 30 dh(4)=1/At*(-Kt(4)*sqrt(h(4)) + qv(4)/3600); 31 32 dz=dh; Código 4.1 Función del modelo no lineal ideal. En este caso, será de vital importancia distinguir entre variables absolutas e incrementales. Para el modelo lineal, no teníamos ningún problema ya que este funciona mediante variables incrementales al igual que los MPC que genera SPCIES. Ahora, el modelo no lineal entiende exclusivamente las variables en forma absoluta, por lo que en TIA PORTAL debemos realizar un ajuste en el bloque principal de esta manera: Figura 4.1 Ajuste para el modelo no lineal en el bloque principal. Como el PLC será el encargado de calcular las acciones de control según el MPC de SPCIES, lo que requiere variables incrementales, tendremos que, en primer lugar, restar el punto de funcionamiento hop al estado siguiente que nos ha enviado MATLAB ® para obtener el incremento del estado deltax con el que trabajará el PLC, y una vez tenemos calculada la acción de control incremental deltau, sumarle a esta el caudal en el punto de operación uop, para enviar a MATLAB ® la acción de control absoluta con la que llamaremos a la función de nuestro modelo no lineal. El bucle de control quedará por tanto así: 4.1 El modelo no lineal simplificado 37 1% Bucle de control MPC 2for k = 1:Nsim 3% Cambios de referencia y cálculo de puntos de equilibrio 4if k >= (Nsim/3) && k < (2*Nsim/3) 5ref = [0.1; 0.2]; 6elseif k >= (2*Nsim/3) 7ref = [-0.1; -0.2]; 8else 9ref = [-0.2; -0.2]; 10 end 11 % Recalcular punto de equilibrio (estado y entrada) 12 nueq = [Ad-eye(4),Bd; Cd, Dd] \ [zeros(4, 1); ref]; 13 xeq = nueq(1:4); 14 ueq = nueq(5:6); 15 16 % Enviamos al MPC 17 for i = 1:4 18 writeValue(opcUaClient, nodosEstadoEq(i), double(xeq(i))); 19 end 20 for i = 1:2 21 writeValue(opcUaClient, nodosAccionEq(i), double(ueq(i))); 22 end 23 24 % Enviar estado actual h al PLC 25 estadocalculado = 1; 26 writeValue(opcUaClient, nodoEstadoCalculado, estadocalculado); 27 for i = 1:4 28 writeValue(opcUaClient, nodosEstado(i), double(h(i))); 29 end 30 pause(0.1); 31 32 % Indicamos al PLC que prosiga con el calculo MPC 33 estadocalculado = 2; 34 writeValue(opcUaClient, nodoEstadoCalculado, estadocalculado); 35 % Esperamos acción calculada del PLC 36 while(readValue(opcUaClient, nodoAccionCalculada)==2) 37 end 38 39 % Leemos accion de control absoluta del PLC 40 if (readValue(opcUaClient, nodoAccionCalculada)==1) 41 for i = 1:2 42 uk(i) = readValue(opcUaClient, nodosAccion(i)); 43 end 44 end 45 46 % Simulación del sistema no lineal 47 modeloIdeal = @(t, z) Modelo4TIdeal(z, uk); 48 [~, z_temp] = ode45(modeloIdeal, [0 Tm], h); 49 h = z_temp(end, :)’; % Nuevo estado 50 51 % Almacenamiento de históricos 52 Mxlqr(:, k + 1) = h; 53 Qx(:, k + 1) = uk; 54 end Código 4.2 Bucle de control para el modelo no lineal. 44 Capítulo 4. Control MPC y simulación de la planta no lineal 4.3.3 Implementación de los controladores PID A continuación, entraremos en el sub bucle mencionado anteriormente, en el que nos encargaremos de que los PID se actualicen a una frecuencia mayor que el MPC, de manera que, leyendo y desescalando los valores recibidos por estos en una apertura de la válvula normalizada entre 0 y 1, integraremos el modelo real de la planta para obtener los caudales actuales, cerrando el lazo de control realimentando los PIDs en la comparación de caudal referencia y caudal actual. 1% Subdivisión del tiempo muestreo Tm, sub bucle 2subTm = 0.25; 3N_sub = Tm / subTm; % 5 / 0.25 = 20 subpasos 4 5for j = 1:N_sub 6 7% Leer caudales del PID (en %) - convertir a m^3/h 8qv = [ 9readValue(opcUaClient, nodoCaudalV1C); 10 readValue(opcUaClient, nodoCaudalV2C); 11 readValue(opcUaClient, nodoCaudalV3C); 12 readValue(opcUaClient, nodoCaudalV4C); 13 ]; 14 qv = (qv / 100) * caudal_max; 15 16 % Calcular apertura normalizada 17 x = qv / caudal_max; 18 x=min(max(x, 0), 1); % Saturamos entre 0 y 1 19 20 % Integramos el modelo con paso subTm 21 modeloReal = @(t, z) Modelo4TReal(z, x); 22 [~, z_temp] = ode45(modeloReal, [0 subTm], z); 23 z = z_temp(end, :)’; 24 25 % Obtener caudales actuales reales 26 [~, qv_actual, ~] = Modelo4TReal(z, x); 27 28 % Escalar caudales actuales reales a porcentaje [0--100] 29 qv_actual_esc = (qv_actual / caudal_max) * 100; 30 31 % Enviar caudales reales al PLC (feedback del PID) 32 writeValue(opcUaClient, nodoCaudalV1A, qv_actual_esc(1)); 33 writeValue(opcUaClient, nodoCaudalV2A, qv_actual_esc(2)); 34 writeValue(opcUaClient, nodoCaudalV3A, qv_actual_esc(3)); 35 writeValue(opcUaClient, nodoCaudalV4A, qv_actual_esc(4)); 36 37 pause(0.2); % tiempo para escribir las variables en el servidor 38 end 39 40 % Extraer las alturas al final del ciclo de control 41 xlqr = z(end-3:end); 42 43 end Código 4.4 Retroalimentación de los PID. El resultado que obtenemos con los parámetros conseguidos mediante el IMC tuning no son nada precisos como podemos apreciar: 4.4 Método alternativo: Sincronización por eventos 45 Figura 4.7 Control de la planta real con parámetros IMC. Para una referencia re f = [−0.2,−0.2]; , se observa un seguimiento de esta con respecto al punto de funcionamiento algo errático, lo cual no ocurría en experimentos anteriores. Además, se da comunmente, en casi todas las pruebas que se realizan, un error de desconexión entre MATLAB ® y el servidor OPC UA. En este caso, poseemos más del doble de variables compartidas en el servidor OPC UA de manera que este se puede ver mucho más afectado ante la saturación de datos de entrada y salida, además necesitando los PID frecuencias mayores. En TIA PORTAL, la tasa de publicación y la tasa de escritura del servidor OPC UA son configurables, pero un valor adecuado es 0.5 segundos. Se han realizado algunas pruebas con el PLC físico en el que no tenemos este problema de desconexión, por lo que debe ser problema o bien del pórtatil empleado, o de el PLCSIM Advanced. 4.4 Método alternativo: Sincronización por eventos Para depender en la menor medida de lo posible de los tiempos de servidor OPC UA, interrupciones cíclicas en el PLC y pausas en MATLAB para intentar conseguir una sincronización temporal, que hemos podido comprobar que se complica, se realizará una sincronización por evento. Para ello, se ha incorporado el código para PIDs discretos en el mismo bloque MPC. Con el fin de respetar los ciclos del PID por cada ciclo del MPC, se ha propuesto lo siguiente: una variable tipo flag denominada EntraMPC, que marcará la entrada en el algoritmo MPC. La inicializaremos a 1, para que en el primer ciclo pueda entrar el algoritmo, y en el inicio del sub bucle de retroalimentación de los PID se pondrá a 0. Una vez terminado este se pondrá a 1, para que pueda volver a entrar el algoritmo MPC. De esta manera, el MPC calculará los caudales referencia una vez por cada número de ciclos especificado para el sub bucle PID. El código PID discreto empleado por ejemplo para la válvula 1 será: // ===== PI VÁLVULA 1 ===== "Data_block_1".e[0] := #qv1 - #qv1A; // Integración: evitar acumulación si el error no es significativo IF "Data_block_1".e[0] > 1.5 OR "Data_block_1".e[0] < -1.5 THEN "Data_block_1".int1 := "Data_block_1".int1 + (#Ts / "Data_block_1".Ti1) * " Data_block_1".e[0]; ELSE "Data_block_1".int1 := "Data_block_1".int1; END_IF; 46 Capítulo 4. Control MPC y simulación de la planta no lineal // Acción de control #qv1C := "Data_block_1".Kp1 * ( "Data_block_1".e[0] + "Data_block_1".int1 ); // Saturar la salida a [0, 100] mas antiwindup IF #qv1C > 100.0 THEN #qv1C := 100.0; "Data_block_1".int1 := "Data_block_1".int1 - (#Ts / "Data_block_1".Ti1) * " Data_block_1".e[0]; ELSIF #qv1C < 0.0 THEN #qv1C := 0.0; "Data_block_1".int1 := "Data_block_1".int1 - (#Ts / "Data_block_1".Ti1) * " Data_block_1".e[0]; END_IF; // Guardar error anterior "Data_block_1".ea[0] := "Data_block_1".e[0]; Código 4.5 Código PID discreto. Siendo el error, el caudal referencia menos el caudal de apertura, la acción de control una acción PI, ya que realmente para las válvulas no hará falta la acción derivativa, e incorporaremos un antiwindup con el fin de evitar oscilaciones debido a la acumulación del término integral. Para el resto de válvulas, simplemente repetiremos este código. El código de MATLAB ® poseerá únicamente la modificación del flag EntraMPC con respecto al comentado en el apartado anterior, por lo que no valdrá la pena comentarlo. En cuanto a TIA PORTAL, será importante crear en el mismo bloque MPC las entradas EntraMPC, qv1 a qv4 para los caudales referencia y qv1A a qv4A para los caudales actuales medidos, y qv1C a qv4C, caudales calculados por los PI, para las salidas. Adjuntaremos un ejemplo de control PID sin antiwindup y otro con antiwindup. Figura 4.8 Control MPC y PIDs sin antiwindup. 4.5 Simulación de la planta no lineal real 47 Figura 4.9 Control MPC y PIDs con antiwindup. Por algún motivo que no se ha llegado a solucionar, el seguimiento con antiwindup no llega a estabilizarse en la referencia, aunque si que reduce las oscilaciones comparado con el PI sin antiwindup, que si que se estabiliza en la referencia deseada ref = [-0.1, -0.1] pero presentando elevadas oscilaciones. 4.5 Simulación de la planta no lineal real Aunque no ha podido ser llevado a cabo el control del modelo no lineal más aproximado a la planta real, la idea de la simulación con WinCC era la siguiente: Figura 4.10 Simulación y entorno HMI de la planta real. En primer lugar como mostramos en la simulación de la planta ideal, la representación en cada instante de muestreo de la altura del agua en cada uno de los tanques con la adición en este caso de la posibilidad de cambio manual de las consignas de altura de los tanques 1 y 2. Por otra parte, la incorporación de las válvulas proporcionales en lugar de las válvulas de tres vías que teníamos en la planta ideal, junto con la referencia o set point (SP) de apertura de la válvula que calcula el MPC y el valor de salida o process value (PV) de los PIDs que deberá tender al SP. Finalmente, un panel 48 Capítulo 4. Control MPC y simulación de la planta no lineal donde poder cambiar los PID de modo automático a modo manual, para introducir manualmente las consignas en los campos SP, y donde probar durante la simulación distintos parámetros de los PID o PI implementados. Con todo esto conseguiríamos el objetivo de la puesta en marcha del sistema por un operario, pudiendo obtener resultados con distintas consignas, parámetros de los PIDs y alternancia entre modo manual y automático. 5 Conclusiones y líneas futuras Quedan fuera del alcance de este trabajo cuestiones como la realización completa y conexionado de un entorno HMI y SCADA completo en WinCC que nos permita controlar la planta de los cuatro tanques a través de nuestro ordenador, pudiendo seleccionar nosotros las referencias que deseamos alcanzar y los parámetros PI o PID para realizar diversos experimentos y optimizar el control de la planta real, ya que se ha planteado exclusivamente el diseño sin una prueba de su funcionamiento, al no haberse logrado el control del modelo no lineal real. También queda pendiente la transición del modo manual al automático y la asignación de diferentes modos como el de emergencia, marcha o paro, siguiendo por ejemplo la guía GEMMA. Además, sería positivo explorar un mayor número de funcionalidades de SPCIES como la implementación de los algoritmos FISTA o EADMM o la posibilidad de poder cambiar las matrices A y B en tiempo real, así como la recogida explícita de distintos datos que ofrece SPCIES como el tiempo de computación o número de iteraciones que se realizan para hallar las soluciones óptimas. Todo esto nos brinda la oportunidad de obtener un conocimiento más profundo de nuestros sistemas. Por último, el objetivo más ambicioso de este proyecto hubiese sido la implementación de los programas desarrollados en el PLC físico y la realización del conexionado, con el propósito de haber puesto por primera vez la planta de los cuatro tanques en funcionamiento controlada con un regulador MPC. 49 Apéndice A Códigos MATLAB® A.1 Código control LQR con observador y cancelación de perturbaciones 1clear 2close all 3%Problema de los 4 tanques para qa0=1.948 y qb0=2 4Tm = 5;%tiempo muestreo 5Tsim =5000;%tiempo total de la simulacion 6Nsim = Tsim/Tm; 7 8 9Mx=zeros(4,Nsim+1); 10 Q = 0.9*eye(4); 11 R = 0.7*eye(2); 12 13 A = 0.03; 14 a1 = 1.3104e-4; 15 a2 = 1.5074e-4; 16 a3 = 9.2673e-5; 17 a4 = 8.8164e-5; 18 gama = 0.3; 19 gamb = 0.4; 20 g=9.81; 21 22 %Punto funcionamiento 23 qa = 1.948; 24 qb = 2; 25 h30 = ((1-gamb)^2)*(qb^2)/((3600^2)*(a3^2)*2*g); 26 h40 = ((1-gama)^2)*(qa^2)/((3600^2)*(a4^2)*2*g); 27 h10 = (1/((a1^2)*2*g))*((a3^2)*2*g*h30+(gama^2)*(qa^2)/(3600^2)+2*a3*sqrt(2*g* h30)*gama*qa/3600); 28 h20 = (1/((a2^2)*2*g))*((a4^2)*2*g*h40+(gamb^2)*(qb^2)/(3600^2)+2*a4*sqrt(2*g* h40)*gamb*qb/3600); 29 u = [qa; qb]; %entrada constante 30 ueq = [qa;qb]; 31 uk = [qa;qb]; 32 33 % Niveles iniciales de los tanques (en metros) 34 h0 = [0.5; 0.5; 0.5; 0.5]; 35 hop = [h10; h20; h30; h40]; 51 52 Capítulo A. Códigos MATLAB® 36 x = h0; 37 Mx(:,1)=h0; 38 Mxlqr(:,1)=hop; 39 ref = [0;0]; 40 nueq = zeros(1,6);%vector que contiene x equilibrio y u equilibrio 41 xeq = h0; 42 43 44 d = [0.2;0.2]; %vector perturbación 45 destk = [0.2;0.2];%vector perturbacion estiamda 46 destk1 = [0.2;0.2];%vector perturbacion estimada anterior 47 48 zk = [x;d]; %Vector de estado extendido que incluye las 4 alturas y perturbacion en la salida 49 zk1 = [x;d]; %6x1 50 51 zestk = [x;d]; %Estimacion del vector de estado 52 yk = x(1:2); 53 54 %Entrada para cancelar el offset: deltau = Kobs*deltax + Kobs*error; 55 56 57 %Matrices A y B en tiempo CONTINUO resultantes de la linealización 58 Ac = [-a1*g/(A*sqrt(2*g*h10)), 0, a3*g/(A*sqrt(2*g*h30)), 0; 59 0, -a2*g/(A*sqrt(2*g*h20)), 0 , a4*g/(A*sqrt(2*g*h40)); 60 0, 0, -a3*g/(A*sqrt(2*g*h30)), 0; 61 0 , 0, 0, -a4*g/(A*sqrt(2*g*h40))]; 62 Bc = [gama/(A*3600), 0; 63 0, gamb/(A*3600); 64 0, (1-gamb)/(A*3600); 65 (1-gama)/(A*3600), 0]; 66 Cc = [1 , 0 , 0, 0; 0, 1, 0, 0]; 67 Dc=zeros(2,2); 68 69 sisc = ss(Ac,Bc,Cc,Dc);%sistema continuo 70 71 %conversion matrices a tiempo DISCRETO mantenedor orden cero 72 sisd = c2d(sisc, Tm); 73 Ad = sisd.A; 74 Bd = sisd.B; 75 Cd = sisd.C; 76 Dd = sisd.D; 77 78 %matrices extendidas para observacion de estados 79 Aext = [Ad , zeros(4,2) ; zeros(2,4) , eye(2)]; 80 Bext = [Bd ; zeros(2,2)]; 81 Cext = [Cd, eye(2)]; 82 83 %calculo de la ganancia por ubicacion de polos 84 polos=[0.5,0.6,0.65,0.8,0.85,0.95]; %garantiza estabilidad al ser menores que 1 85 86 87 %Hallamos la ganancia de observabilidad 88 L = place(Aext’, Cext’, polos)’; %con comando place podemos hallar tanto ganancia de 89 %observabilidad como de controlabilidad, como son opuestas pillamos la A.1 Código control LQR con observador y cancelación de perturbaciones 53 90 %traspuesta de Aext y Cext y trasponemos la ganancia para dimensiones OK 91 92 tiempo = 0:Tsim; 93 94 %%Control LQR 95 96 %matriz de coste del estado Q y matriz coste entrada R 97 [K]=dlqr(Ad,Bd,Q,R); 98 99 %xeq=[Ad-eye(4),B; Cd, Dd]\[zeros(4, 1);ref]; 100 101 d=[0.2;0.2]; 102 103 Ref=zeros(2,1); 104 xlqr=zeros(4,1); 105 zestk1=zeros(6,1); 106 uk1=zeros(2,1); 107 yk1=zeros(2,1); 108 109 110 % Bucle de control LQR con perturbacion y observador 111 for k = 1:Tsim 112 113 if k >= (Tsim/3) && k < (2*Tsim/3) 114 ref = [1; 0.9]; 115 116 elseif k >= (2*Tsim/3) 117 ref = [0.7; 0.6]; 118 119 120 else 121 ref = [0.5; 0.5]; % Punto de operación 122 end 123 if k == (Tsim/2) 124 d = 3* [0.1;0.1]; 125 end 126 127 % Fórmula LQR, a partir de Tsim/3 entra en funcionamiento el control 128 if k >= (Tsim/3) 129 zestk = Aext * zestk1 + Bext * uk1 + L * (yk1 - Cext * zestk1); 130 xestk=zestk(1:4); 131 destk = zestk(5:6); 132 133 nueq = [Ad-eye(4),Bd; Cd, Dd]\[zeros(4, 1);ref - destk];%calculamos alturas equilibrio y caudales equilibrio para esta referencia 134 135 xd = nueq(1:4); 136 ud = nueq(5:6); 137 138 uk = ud - K*(xestk - xd); 139 140 % Almacenamiento de los Estados pasados 141 % Antes de la predicción 142 x1 = xlqr; 143 uk1 = uk; 144 yk1 = yk; 60 Capítulo B. Códigos TIA PORTAL #lambda[#ii, #jj] := 0.0; END_FOR; END_FOR; #k_iter := 0; // Reseteo del contador de iteraciones al recibir nuevo estado // Actualizar los primeros #n elementos de #b FOR #jj := 0 TO #n-1DO #b[#jj] := 0.0; // Inicializa #b[#j] a 0 FOR #ii := 0 TO #n-1DO #b[#jj] := #b[#jj] - #AB[#jj, #ii] * #x0[#ii]; // Calcula #b[#j] END_FOR; END_FOR; // Actualizar la referencia FOR #jj := 0 TO #n-1DO #q_1[#jj] := #Q[#jj] * #xr[#jj]; // Calcula #q[#j] #qT[#jj] := 0.0; // Inicializa #qT[#j] a 0 FOR #ii := 0 TO #n-1DO #qT[#jj] := #qT[#jj] + #T[#jj, #ii] * #xr[#ii]; // Calcula #qT[#j] END_FOR; END_FOR; FOR #jj := 0 TO #m-1DO #q_1[#jj + #n] := #R[#jj] * #ur[#jj]; // Calcula #q[#j + #nn_] END_FOR; END_IF; IF #new_data_available = 2 AND #new_state_received = 2 THEN //si recibimos un nuevo estado empezamos optimizacion FOR #iter_count := 0 TO #iter_per_cycle - 1 DO #k_iter := #k_iter + 1; // Paso 0: Guardar valor de v en la variable v1 // memcpy(v1_0, v_0, sizeof(double)*mm_); //Copiamos un bloque de memoria de una ubicacion a otra, en este caso primer elemento de v FOR #ii := 0 TO #m-1DO #v1_0[#ii] := #v_0[#ii]; END_FOR; // memcpy(v1, v, sizeof(double)*(NN_-1)*nm_); //copiamos elementos de v en otro vector FOR #ii := 0 TO #NN-2DO FOR #jj := 0 TO #nm-1DO B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 61 #v1[#ii, #jj] := #v[#ii, #jj]; END_FOR; END_FOR; // memcpy(v1_N, v_N, sizeof(double)*nn_); //copiamos elemento final de v en otro vector FOR #ii := 0 TO #n-1DO #v1_N[#ii] := #v_N[#ii]; END_FOR; // Step 1: Minimizar con respecto a z // Calculo del vector q_hat = lambda - rho*v //Guardamos vector q_hat en z para ahorrar memoria y computacion // Calculo de los primeros m elementos FOR #jj := 0 TO #m-1DO #z_0[#jj] := #q_1[#jj + #n] + #lambda_0[#jj] - #rho * #v_0[#jj]; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 0 TO #NN-2DO FOR #jj := 0 TO #nm-1DO #z[#ll, #jj] := #q_1[#jj] + #lambda[#ll, #jj] - #rho * #v[#ll, # jj]; END_FOR; END_FOR; // Calculo de los elementos finales FOR #jj := 0 TO #n-1DO #z_N[#jj] := #qT[#jj] + #lambda_N[#jj] - #rho * #v_N[#jj]; END_FOR; //Calculo del lado derecho del sistema de ecuaciones Wc: -G’*H_hat^(-1)* q_hat - b // La guardamos en mu para ahorrar algo de memoria // Calculo de los primeros n elementos FOR #jj := 0 TO #n-1DO #mu[0, #jj] := #Hi[0, #jj] * #z[0, #jj] - #b[#jj]; FOR #ii := 0 TO #m-1DO 62 Capítulo B. Códigos TIA PORTAL #mu[0, #jj] := #mu[0, #jj] - #AB[#jj, #ii + #n] * #Hi_0[#ii] * # z_0[#ii]; END_FOR; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 1 TO #NN-2DO FOR #jj := 0 TO #n-1DO #mu[#ll, #jj] := #Hi[#ll, #jj] * #z[#ll, #jj]; FOR #ii := 0 TO #nm-1DO #mu[#ll, #jj] := #mu[#ll, #jj] - #AB[#jj, #ii] * #Hi[#ll - 1, #ii] * #z[#ll - 1, #ii]; END_FOR; END_FOR; END_FOR; // Calculo de los ultimos n elementos FOR #jj := 0 TO #n-1DO #mu[#NN - 1, #jj] := 0.0; FOR #ii := 0 TO #n-1DO #mu[#NN - 1, #jj] := #mu[#NN - 1, #jj] + #Hi_N[#jj, #ii] * #z_N [#ii]; END_FOR; FOR #ii := 0 TO #nm-1DO #mu[#NN - 1, #jj] := #mu[#NN - 1, #jj] - #AB[#jj, #ii] * #Hi[#NN - 2, #ii] * #z[#NN - 2, #ii]; END_FOR; END_FOR; // Calculo de mu, solucion del sistema de ecuaciones W*mu = -G’*H^(-1)* q_hat - beq // FORWARD SUBSTITUTION FOR #jj := 0 TO #n-1DO FOR #ii := 0 TO #jj-1DO #mu[0, #jj] := #mu[0, #jj] - #Beta[0, #ii, #jj] * #mu[0, #ii]; B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 63 END_FOR; #mu[0, #jj] := #Beta[0, #jj, #jj] * #mu[0, #jj]; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 1 TO #NN-2DO FOR #jj := 0 TO #n-1DO FOR #ii := 0 TO #n-1DO #mu[#ll, #jj] := #mu[#ll, #jj] - #Alpha[#ll - 1, #ii, #jj] * #mu[#ll - 1, #ii]; END_FOR; FOR #ii := 0 TO #jj-1DO #mu[#ll, #jj] := #mu[#ll, #jj] - #Beta[#ll, #ii, #jj] * #mu[# ll, #ii]; END_FOR; #mu[#ll, #jj] := #Beta[#ll, #jj, #jj] * #mu[#ll, #jj]; END_FOR; END_FOR; // Calculo de los ultimos n elementos FOR #jj := 0 TO #n-1DO FOR #ii := 0 TO #n-1DO #mu[#NN - 1, #jj] := #mu[#NN - 1, #jj] - #Alpha[#NN - 2, #ii, # jj] * #mu[#NN - 2, #ii]; END_FOR; FOR #ii := 0 TO #jj-1DO #mu[#NN - 1, #jj] := #mu[#NN - 1, #jj] - #Beta[#NN - 1, #ii, #jj ] * #mu[#NN - 1, #ii]; END_FOR; #mu[#NN - 1, #jj] := #Beta[#NN - 1, #jj, #jj] * #mu[#NN - 1, #jj]; END_FOR; // BACKWARD SUBSTITUTION // Calculo de los ultimos nn_ elementos 64 Capítulo B. Códigos TIA PORTAL // Calculo de los ultimos n elementos FOR #jj := #n - 1 TO 0 BY -1 DO FOR #ii := #n - 1 TO #jj + 1 BY -1 DO #mu[#NN - 1, #jj] := #mu[#NN - 1, #jj] - #Beta[#NN - 1, #jj, #ii ] * #mu[#NN - 1, #ii]; END_FOR; #mu[#NN - 1, #jj] := #Beta[#NN - 1, #jj, #jj] * #mu[#NN - 1, #jj]; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := #NN - 2 TO 1 BY -1 DO FOR #jj := #n - 1 TO 0 BY -1 DO FOR #ii := #n - 1 TO 0 BY -1 DO #mu[#ll, #jj] := #mu[#ll, #jj] - #Alpha[#ll, #jj, #ii] * #mu [#ll + 1, #ii]; END_FOR; FOR #ii := #n - 1 TO #jj + 1 BY -1 DO #mu[#ll, #jj] := #mu[#ll, #jj] - #Beta[#ll, #jj, #ii] * #mu[# ll, #ii]; END_FOR; #mu[#ll, #jj] := #Beta[#ll, #jj, #jj] * #mu[#ll, #jj]; END_FOR; END_FOR; // Calculo de los primeros n elementos FOR #jj := #n - 1 TO 0 BY -1 DO FOR #ii := #n - 1 TO 0 BY -1 DO #mu[0, #jj] := #mu[0, #jj] - #Alpha[0, #jj, #ii] * #mu[1, #ii]; END_FOR; FOR #ii := #n - 1 TO #jj + 1 BY -1 DO #mu[0, #jj] := #mu[0, #jj] - #Beta[0, #jj, #ii] * #mu[0, #ii]; END_FOR; #mu[0, #jj] := #Beta[0, #jj, #jj] * #mu[0, #jj]; B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 65 END_FOR; // Calcular z (importante, de antes en este punto tendremos z = q_hat) // Calculo de los primeros m elementos FOR #jj := 0 TO #m-1DO FOR #ii := 0 TO #n-1DO #z_0[#jj] := #z_0[#jj] + #AB[#ii, #jj + #n] * #mu[0, #ii]; END_FOR; #z_0[#jj] := - #Hi_0[#jj] * #z_0[#jj]; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 0 TO #NN-2DO FOR #jj := 0 TO #n-1DO #z[#ll, #jj] := #z[#ll, #jj] - #mu[#ll, #jj]; END_FOR; FOR #jj := 0 TO #nm-1DO FOR #ii := 0 TO #n-1DO #z[#ll, #jj] := #z[#ll, #jj] + #AB[#ii, #jj] * #mu[#ll + 1, # ii]; END_FOR; #z[#ll, #jj] := - #Hi[#ll, #jj] * #z[#ll, #jj]; END_FOR; END_FOR; // Calculo de los ultimos n elementos FOR #jj := 0 TO #n-1DO #aux_N[#jj] := #z_N[#jj] - #mu[#NN - 1, #jj]; END_FOR; FOR #jj := 0 TO #n-1DO #z_N[#jj] := 0.0; FOR #ii := 0 TO #n-1DO 66 Capítulo B. Códigos TIA PORTAL #z_N[#jj] := #z_N[#jj] - #Hi_N[#jj, #ii] * #aux_N[#ii]; END_FOR; END_FOR; // Paso 2: Minimizar con respecto a v // Calculo de las primeras mm variables FOR #jj := 0 TO #m-1DO #v_0[#jj] := #z_0[#jj] + #rho_i * #lambda_0[#jj]; //aplicamos restricciones IF #v_0[#jj] < #LB[#jj + #n] THEN #v_0[#jj] := #LB[#jj + #n]; END_IF; IF #v_0[#jj] > #UB[#jj + #n] THEN #v_0[#jj] := #UB[#jj + #n]; END_IF; END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 0 TO #NN-2DO FOR #jj := 0 TO #nm-1DO #v[#ll, #jj] := #z[#ll, #jj] + #rho_i * #lambda[#ll, #jj]; IF #v[#ll, #jj] <= #LB[#jj] THEN #v[#ll, #jj] := #LB[#jj]; ELSIF #v[#ll, #jj] >= #UB[#jj] THEN #v[#ll, #jj] := #UB[#jj]; END_IF; END_FOR; END_FOR; // Calculo de los ultimos n elementos FOR #jj := 0 TO #n-1DO #v_N[#jj] := #z_N[#jj] + #rho_i * #lambda_N[#jj]; B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 67 IF #v_N[#jj] > #UB[#jj] THEN #v_N[#jj] := #UB[#jj]; ELSIF #v_N[#jj] < #LB[#jj] THEN #v_N[#jj] := #LB[#jj]; END_IF; END_FOR; // Paso 3: Actualizar lambda // Calculo de los primeros mm elementos FOR #jj := 0 TO #m-1DO #lambda_0[#jj] := #lambda_0[#jj] + #rho * (#z_0[#jj] - #v_0[#jj]); END_FOR; // Calculo del resto de elementos excepto los n ultimos FOR #ll := 0 TO #NN-2DO FOR #jj := 0 TO #nm-1DO #lambda[#ll, #jj] := #lambda[#ll, #jj] + #rho * (#z[#ll, #jj] - #v[#ll, #jj]); END_FOR; END_FOR; // Calculamos ultimos n elementos FOR #jj := 0 TO #n-1DO #lambda_N[#jj] := #lambda_N[#jj] + #rho * (#z_N[#jj] - #v_N[#jj]); END_FOR; END_FOR; // Paso 4: Calculo del residuo #res_flag := 0; // Reset del flag residuo // Calcular primeros m elementos FOR #jj := 0 TO #m-1DO #res_fixed_point := #v1_0[#jj] - #v_0[#jj]; #res_primal_feas := #z_0[#jj] - #v_0[#jj]; // Obtener valores absolutos 68 Capítulo B. Códigos TIA PORTAL IF #res_fixed_point < 0.0 THEN #res_fixed_point := - #res_fixed_point; END_IF; IF #res_primal_feas < 0.0 THEN #res_primal_feas := - #res_primal_feas; END_IF;// Fin obtencion valores absolutos IF #res_fixed_point > #tol_d OR #res_primal_feas > #tol_p THEN #res_flag := 1; EXIT; END_IF; END_FOR; // Calculamos ultimos n elementos IF #res_flag = 0 THEN FOR #jj := 0 TO #n-1DO #res_fixed_point := #v1_N[#jj] - #v_N[#jj]; #res_primal_feas := #z_N[#jj] - #v_N[#jj]; // Obtencion valores absolutos IF #res_fixed_point < 0.0 THEN #res_fixed_point := - #res_fixed_point; END_IF; IF #res_primal_feas < 0.0 THEN #res_primal_feas := - #res_primal_feas; END_IF;// Fin obtencion valores absolutos IF #res_fixed_point > #tol_d OR #res_primal_feas > #tol_p THEN #res_flag := 1; EXIT; END_IF; END_FOR; END_IF; // Calculamos el resto de elementos IF #res_flag = 0 THEN FOR #ll := 0 TO #NN-2DO B.1 Código MPC con algoritmo ADMM traducido a SCL para TIA PORTAL 69 FOR #jj := 0 TO #nm-1DO #res_fixed_point := #v1[#ll, #jj] - #v[#ll, #jj]; #res_primal_feas := #z[#ll, #jj] - #v[#ll, #jj]; // Obtener valores absolutos IF #res_fixed_point < 0.0 THEN #res_fixed_point := - #res_fixed_point; END_IF; IF #res_primal_feas < 0.0 THEN #res_primal_feas := - #res_primal_feas; END_IF;// Fin obtencion valores absolutos IF #res_fixed_point > #tol_d OR #res_primal_feas > #tol_p THEN #res_flag := 1; EXIT; END_IF; END_FOR; IF #res_flag = 1 THEN EXIT; END_IF; END_FOR; END_IF; // Paso 5: Condicion de salida IF #res_flag = 0 THEN #new_data_available := 1; #e_flag := 1; ELSIF #k_iter >= #k_max THEN #new_data_available := 1; #e_flag := -1; END_IF; END_IF; IF #new_data_available = 1 THEN #k_out := #k_iter; Glosario 77 Glosario ETSI Escuela Técnica Superior de Ingeniería. 77 HMI Interfaz Hombre-Máquina. 77 LQR Regulador Lineal Cuadrático. 77 MPC Control Predictivo Basado en Modelo. 77 OPC UA Open Platform Communications Unified Architecture. 77 PID Proporcional-Integral-Derivativo. 77 PLC Controlador Lógico Programable. 77 SCADA Supervisión, Control y Adquisición de Datos. 77 TIA Portal Totally Integrated Automation Portal. 77 US Universidad de Sevilla. 77 WinCC Windows Control Center. 77 79