Desarrollo de solvers para MPC en el lenguaje de programación Julia

Equation Chapter 1 Section 1 Trabajo Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Desarrollo de solvers para MPC en el lenguaje de programación Julia Autor: Rubén León Fuentes Tutores: Pablo Krupa García Ignacio Alvarado Aldea Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 iii Trabajo Fin de Grado Ingeniería Electrónica, Robótica y Mecatrónica Desarrollo de solvers para MPC en el lenguaje de programación Julia Autor: Rubén León Fuentes Tutores: Pablo Krupa García Ignacio Alvarado Aldea Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Trabajo Fin de Grado: Desarrollo de solvers para MPC en el lenguaje de programación Julia Autor: Rubén León Fuentes Tutores: Pablo Krupa García Ignacio Alvarado Aldea El tribunal nombrado para juzgar el proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: El Secretario del Tribunal Sevilla, 2021 vii AGRADECIMIENTOS A mis compañeros, los que estaban al principio y los que encontré por el camino. A Marina, por el apoyo y la alegría que me ayudan a avanzar. A mi familia, que me ha traído hasta aquí. A Ana, por ser una fuente inmensa de cariño y compañía durante tantos y tantos días. A Pablo, sin cuya inestimable ayuda este proyecto no sería ni la mitad de lo que es. Todo lo que jamás se ha escrito lo ha escrito un único autor. RESUMEN El control predictivo por modelo (MPC) representa una de las técnicas avanzadas de control que vienen abriéndose paso en el sector desde hace unas décadas. En este proyecto, se desarrolla un solver para la resolución del problema matemático que surge de la implementación de MPC. El solver se desarrolla en el lenguaje novedoso y en alza Julia y utilizada un método de optimización bien conocido en la disciplina como es ADMM. Una vez se ha expuesto el marco teórico y se ha adecuado la formulación al problema que se desea resolver, se desarrolla una implementación del mismo tanto en Matlab como en Julia. Posteriormente, se valida el funcionamiento del solver a través de una serie de simulaciones frente a diferentes tipos de sistemas o parámetros de configuración. ix ABSTRACT Since last decades, some advanced control techniques are breaking through the control systems field. Model predictive control (MPC) represents one of them. In this project, a solver for the mathematical problem arising from MPC implementation is developed. The solver is written in Julia, a new language currently on the rise, and it involves the use of the optimization method ADMM. After discussing the theoretical framework and adjusting the formulation to the problem at hand, the solver implementation is developed both in Matlab and Julia. Afterwards, the solver performance is validated through a series of simulations with different settings for the controlled systems and system parameters. 1 1 INTRODUCCIÓN El control predictivo basado en modelo (MPC) es un conjunto de técnicas de control basadas en el uso de un modelo para realizar predicciones del estado futuro del sistema a controlar como parte del proceso de control ([1], [2]). Supone el marco en el que se desarrolla este proyecto y es el punto central sobre el que se mueven todas las demás ideas planteadas. Como comienzo de este proyecto se introduce al lector en qué situación se encuentra actualmente la tecnología de MPC en el sector de los sistemas de control. De aquí se deducen las motivaciones y objetivos que se presuponen como guía de las tareas realizadas a lo largo del proyecto. También se incluye un pequeño resumen por capítulos de este documento. 1.1 Estado del arte La aplicación de sistemas de control orientados a la predicción con modelo se comienza a dar en la década de 1970, aunque algunos de los conceptos fundamentales de MPC surgen con anterioridad aplicados a otras técnicas de control. Estos sistemas se abren paso en la industria petroquímica [3]. Utilizan modelos de respuesta a impulso y salto escalonado y ya incluyen la implementación de restricciones en sistemas multivariables, sin embargo, no cuentan con una teoría desarrollada que los sustente y los provea de resultados en cuanto a estabilidad o robustez [1]. De manera independiente surgen controladores en basados modelos de función de transferencia en sistemas monovariable que aplican ideas de control adaptativo. Son varios los diseños de controlador que aparecen bajo esta premisa [1]. También aparecen bajo este contexto controladores basados en modelo de espacio de estados, que permiten fácilmente ampliar a sistemas multivariable o analizar la estabilidad [4]. Aun así, la teoría no alcanzaba a proponer resultados generales de estabilidad para los controladores que utilizaban horizonte de predicción finito. En la década de los noventa, con la aparición de nuevos controladores centrados en espacio de estados y enfocados en este asunto, se le puso solución [5]. Actualmente, las nuevas tecnologías que están surgiendo en el contexto del MPC se centran en obtener controladores robustos. Entre la comunidad, hoy en día, las técnicas de control que se engloban bajo MPC son un estándar aceptado para sistemas lineales y relativamente lentos. Los sistemas no lineales o sometidos a tiempos de muestreo muy bajos representan el futuro del sector. A pesar de la gran cantidad de controladores desarrollados e investigados en el mundo académico y del interés que producen desde el punto de vista teórico, las ideas comentadas en el último párrafo han frenado históricamente la penetración de estas tecnologías en la industria. Entre los sectores dónde más recurrente es su uso destaca el ya mencionado como origen, sector petroquímico [6]. También se puede encontrar en la industria aeroespacial o automovilística. Sectores donde la no linealidad de los sistemas se encuentra más presente suponen un reto para esta técnica de control, aunque se han realizado progresos en los últimos años. Este proyecto se enfoca en particular en un hito relativamente reciente en la historia de MPC: su aplicación en sistemas embebidos. Esto representa un paradigma útil que permitiría llevar técnicas avanzadas de control y sus ventajas al tipo de dispositivos controladores más presentes en la industria, como son, por ejemplo, los Programmable Logic Controllers (PLC). En esta línea de investigación destaca la aparición de sistemas lineales y no robustos, aunque también aparecen casos contrarios [7]. En ella se encuentran dos formas diferentes de aplicar MPC: MPC online y MPC explícito. El primero se corresponde con las técnicas de control predictivo en las que el problema de optimización que surge de MPC se resuelve en cada tiempo de muestreo del sistema [8]. Existen solvers diseñados específicamente con este propóstio en mente, lo que permite resultados muy eficientes. Algunos ejemplos son μAO-MPC [9], HPMPC [10] o SPCIES [11]. La otra forma de aplicarlo, la explícita, se aprovecha de que la acción de control óptima de un MPC lineal estándar responde a un mapa afín a trozos [12]. Introducción 2 Este mapa se puede calcular en un dispositivo ajeno al controlador, dando lugar a una look-up table. Esta tabla permite obtener la acción de control sin resolver el problema de optimización. Esta variante destaca por su velocidad, sin embargo, las dimensiones de la tabla crecen rápidamente con las dimensiones del problema. Esto provoca que se requiera mucha memoria y tiempo para encontrar la solución, tanto que solo es aplicable a sistemas pequeños con pocas restricciones. 1.2 Motivación del proyecto El control predictivo basado en modelo es una técnica relativamente avanzada y eficiente de control que, a pesar de estar ampliamente estudiada en el ámbito académico, no tiene toda la penetración que cabría esperar en la industria. Uno de los motivos a resaltar que provocan esta situación es la alta capacidad computacional que requieren los dispositivos controladores utilizados para implementar esta técnica. Esto se debe a que MPC (online) requiere resolver un problema de optimización matemática en cada tiempo de muestreo. Este factor es crítico en aquellos sistemas con tiempos de muestreo muy bajos, pues los tiempos de computo hacen de cuello de botella para cuan bajo pueden tomarse los tiempos de muestreo. Este proyecto puede enmarcarse en el trabajo desarrollado por el profesor de la Universidad de Sevilla Pablo Krupa, quien tiene a sus espaldas numerosas investigaciones centradas en esto. Krupa, que además es uno de lo supervisores de este proyecto, ha escrito varios artículos académicos al respecto de la implementación de MPC en sistemas de escasos recursos [13]. En ellos se centra en aplicar enfoques más eficientes y competitivos al desarrollo de MPC. Además, ha desarrollado la herramienta SPCIES [11], que permite generar de forma automatizada solvers que implementan MPC en lenguajes como C o Maltab. Siendo el lenguaje de programación donde se implementa el solver una característica importante, se abre la oportunidad de aplicar estos conocimientos a otros lenguajes que permitan un desarrollo aún mayor de la eficiencia. En los últimos años se está asistiendo a una abundancia en la creación de nuevos lenguajes de carácter más específico, con la intención de enfocarlos a tareas más concretas y con cualidades especiales que explotar en determinados campos. Así nace Julia [14], un lenguaje caracterizado por sus logros en lo que a velocidad de ejecución respecta. Ya existen solvers en Julia de propósito general para problemas de optimización convexa, como es el caso de COSMO [15]. Estos solvers pueden adaptarse a MPC, pero no surgen con este en mente. Es así que este proyecto pretende sumar un granito de arena a la implementación de solvers eficientes para MPC a través de poner en juego estas ideas en un programa desarrollado en Julia. 1.3 Objetivos del proyecto El objetivo central es obtener un solver funcional en el lenguaje de programación Julia. El funcionamiento del solver se debe poder validar a través de simulación. El solver debe ser capaz de ejecutar el algoritmo ADMM (Alternating direction method of multipliers) [16] para encontrar la solución al problema quadratic programming (QP) al que se enfrenta, a raíz del desarrollo de la formulación elegida para MPC, en cada tiempo de muesteo. Se pretende que el solver tenga un funcionamiento adecuado frente a diferentes tipos de sistemas. Esto incluye sistemas con diferente número de estados o señales de control. También debe probarse el funcionamiento con diferentes configuraciones de los parámetros que marcan las diferentes tecnologías en juego, como son MPC o ADMM. De forma previa, ese mismo solver se validará en Matlab. Además de comprobar que el solver puede alcanzar la referencia en una simulación, se comparará su comportamiento con el de otros solvers disponibles, como por ejemplo el solver cuadrático nativo de Matlab o COSMO. Se espera una respuesta similar en cuanto a estados y señales de control, pero se considera un reto particularmente difícil alcanzar la calidad de estos solvers. Aun así, se celebrará todo lo que pueda alcanzarse en cuestión de velocidad. Otra cuestión que tratar es la eficiencia. Se desea que el solver tenga un diseño modular que permita reducir al mínimo el código necesario a ejecutar en cada tiempo de muestreo. Se aplica álgebra sparsa en diversas operaciones para reducir los costes de memoria de los datos utilizados y para reducir el número de operaciones a realizar. La implementación de estas cuestiones debe ser validada de modo que el solver eficiente que las aplica tenga un funcionamiento correcto igual al de un solver que no las aplique. De la implementación en Julia se espera que ofrezca unos tiempos competitivos, mejorando los resultados de Matlab y acercándose todo lo posible a los de COSMO. 3 3 Estos representan los objetivos desde un punto de vista práctico. Aun así, no debe olvidarse que este proyecto es un trabajo académico formativo, por lo que también se valoran los conocimientos obtenidos en las diferentes disciplinas que lo conforman. 1.4 Resumen por capítulos El presente documento se encuentra estructurado según se muestra a continuación. Presentados de forma introductoria el estado del arte y la ubicación en el mismo de este proyecto, se continúa explicando las nociones básicas de los principales conceptos que sostienen la teoría que permea este proyecto, todo ello en el Capítulo 2. En el Capítulo 3 se aplican dichos conceptos y se desarrollan, conectándose unos a otros, de forma que se otorga al lector el cuerpo teórico del solver desarrollado. Se incluyen todas las ecuaciones que conforman el proceso de resolución del problema del MPC. El Capítulo 4 consiste en la exposición de la estructura y el diseño del solver, tanto así como las decisiones que florecieron de los retos que supuso la elaboración del mismo. Los resultados obtenidos de este proceso son mostrados al lector a través de los pseudocódigos de cada una de las funciones que componen la herramienta desarrollada. Los resultados de las simulaciones con las que se ha validado el correcto funcionamiento del solver se encuentran en el Capítulo 5. En él están expuestas las gráficas relativas al control de estados y señales de control. También se ilustran las iteraciones y los tiempos de ejecución obtenidos. Los diferentes apartados representan a cada una de las simulaciones, con una miríada de configuraciones diferentes, para que pueda comprobarse el funcionamiento del solver ante situaciones interesantes de la práctica e inducir de ahí las diferentes características de este. Finalmente se incluye una tabla con estadísticas de los tiempos de ejecución para cada una de las simulaciones. Las conclusiones obtenidas de dichas simulaciones son presentadas en el Capítulo 6. Ahí se reflexiona sobre lo que ha supuesto desarrollar este proyecto, además de los beneficios obtenidos de la elaboración del mismo. Se analizan los resultados del capítulo anterior y se desarrollan los éxitos y fracasos que suponen los resultados. 5 2 MARCO TEÓRICO En este apartado se explican los conceptos teóricos sobre los que se fundamenta este proyecto, de la forma más general posible. Es decir, en lo que concierne estas explicaciones, cada uno de los conceptos es independiente de los demás. Por un lado, se comienza por plantear el problema, de lo que se parte y cuáles son los detalles que se deben afrontar. Seguidamente, se habla de la técnica de control que se ha utilizado, como surge esta técnica, por qué se caracteriza, cuál es su alcance y qué ventajas e inconvenientes aporta. La acción de control es calculada por un solver. El algoritmo que implementa el solver es otra de las bases del proyecto. También se responde a en qué consiste el algoritmo, cuáles son sus pasos, qué condiciones debe cumplir, en qué conceptos matemáticos se basa. Por último, se introduce una de las herramientas utilizadas para mejorar la eficiencia del mismo. 2.1 Planteamiento del problema Los sistemas de control están, a grandes rasgos, compuestos de dos elementos principales: el sistema físico que se desea controlar y el controlador encargado de realizar dicha tarea. El grueso del esfuerzo lo ocupa diseñar el controlador, sin embargo, es necesario definir previamente que condiciones son las que describen el sistema con el que se está trabajando. Todo sistema físico es gobernado por unas leyes que en mayor o menor grado pueden expresarse matemáticamente a través de un modelo. Para este proyecto, a lo largo del desarrollo teórico y diseño del sistema, no se ha trabajado con un sistema real concreto, sino sobre un modelo genérico bajo el que pueden describirse diversos sistemas con una serie de características comunes. El modelo al que se ha decidido recurrir es una descripción del sistema en espacio de estados dada por el siguiente modelo lineal e invariante en el tiempo 𝑥(𝑘+1)=𝐴𝑥(𝑘)+𝐵𝑢(𝑘), 𝑥∈ℝ𝑛,𝑢∈ℝ𝑚⁡ (1) donde 𝑥(𝑘) y 𝑢(𝑘) son los estados y señales de control, respectivamente, para el instante 𝑘. Dependiendo del sistema físico que represente pueden usarse unos u otros métodos para obtener el modelo, aunque es común obtenerlo a través de la linealización de las ecuaciones que gobiernan el sistema aplicadas a un determinado punto de operación. El modelo en la descripción en espacio de estados está formado por un par de ecuaciones, sin embargo, como a lo largo de este proyecto solo se trabaja con el control de los estados y señales de control, se obvia la ecuación correspondiente a la salida del sistema. El modelo utilizado incluye unas restricciones, que pueden representar, por ejemplo, límites físicos que puede alcanzar el sistema o zonas de operación del sistema que no se desea que sean sobrepasadas por alguna determinada cuestión. Estas restricciones, tanto para los estados como para las señales de control, se modelan a través de 𝑥≤𝑥≤𝑥, 𝑢≤𝑢≤𝑢 (2) El objetivo de control del problema planteado es llevar el sistema a una referencia constante dada, asumiendo que es un punto de operación del modelo, sin violar las restricciones. No violar las restricciones aplica tanto al proceso del control (en ningún instante del control ni los estados ni las señales de control pueden salir del marco que le imponen las restricciones) como a la referencia que se busca alcanzar (es decir, la referencia también debe cumplir las restricciones). Esto último es lógico, pues sino la referencia sería inalcanzable. Además de las ecuaciones mostradas, a lo largo de los desarrollos realizados en este proyecto se asumen una serie de proposiciones de cara a definir más estrictamente las propiedades del modelo del sistema. i. Es sistema (1) es controlable Marco teórico 6 ii. Las retricciones (2) tienen un interior no vacío, es decir, la restricción inferior debe ser menos que la restricción superior. iii. La referencia (𝑥𝑟,𝑢𝑟) es un punto de equilibrio del sistema, es decir, cumple con la ecuación 𝑥𝑟=𝐴𝑥𝑟+𝐵𝑢𝑟 (3) iv. La referencia (𝑥𝑟,𝑢𝑟) cumple las restricciones (2), es decir, se tiene que 𝑥≤𝑥𝑟≤𝑥, 𝑢≤𝑢𝑟≤𝑢 (4) 2.2 Control predictivo basado en modelo El control predictivo basado en modelo (MPC) es un conjunto de técnicas de control basadas en el uso de un modelo para obtener la acción de control a través de la resolución de un problema de optimización matemática ([1], [2]). El modelo permite realizar predicciones del comportamiento del sistema en instantes futuros. El problema de optimización consiste en la minimización de una función objetivo llamada función de coste. Además de estas dos características, se incluye una tercera: el deslizamiento del horizonte de predicción. Esto significa que, para cada instante en el que se realiza la predicción, se calculan predicciones para un número 𝑁 de instantes futuros, siendo 𝑁 conocido como horizonte de predicción. Para el siguiente instante, se vuelven a culcular 𝑁 predicciones, por lo que se desplaza en una unidad el último instante calculado en la predicción con respecto a la predicción anterior. En el contexto de este proyecto, en el que se trabaja con sistemas discretos, los instantes de tiempo son los intervalos de muestreo. En las predicciones se incluye el cálculo de la acción de control para dichos instantes. En cada tiempo de muestreo, la acción de control predicha para el instante actual (la cual es solución del problema que se plantea resolver) es aplicada al sistema que se desea controlar. Por tanto, lo más importante para las técnicas que siguen la metodología de MPC son el modelo del sistema y la resolución del problema de optimización. La diferencia entre las distintas formulaciones y algoritmos de MPC radica en la elección del tipo de modelo o de la función de coste, principalmente. Este conjunto de características comunes entre todos los diferentes controladores que puede agruparse en la familia de MPC hacen que la técnica de control seguida por estos controladores pueda describirse en los siguientes pasos. En primer lugar, para cada instante de tiempo, se calculan las salidas y estados para tantos instantes como indique el horizonte de predicción. Estos cálculos dependen de dos factores, los estados pasados del sistema y las señales de control futuras, que son aquellas que se desean calcular para posteriormente ser enviadas al sistema. En segunda instancia, las señales de control futuras son calculadas a través de un problema de optimización, que, por lo general, tiene como función objetivo una función cuadrática relacionada con la diferencia entre los estados y señales de control con la referencia en el instante actual. Aunque para casos más simples puede encontrarse una solución analítica, para los casos con restricciones, función objetivo cuadrática y modelo linear, que suponen los más interesantes, pues la adicción de restricciones es un gran potencial de MPC, debe encontrarse la solución a través de métodos numéricos. Finalmente, la señal de control calculada para el instante actual (el instante inicial en las predicciones) es tomada como acción de control y es enviada al sistema. En el nuevo instante, en vez de tomar el calculo que se hizo del estado para ese instante, se lee el estado actual del sistema, que siempre será más fiable; y se repite el proceso desde el primer paso [1]. Para obtener un modelo de predicción pueden utilizarse diversos métodos, ya sean heurísticos o analíticos. Para este proyecto se toma un modelo en espacio de estados. Su principal ventaja es que permite modelar sistemas multivariable de manera sencilla. Se supone que el estado del sistema es accesible y no es necesario un observador [1]. El concepto de utilizar un modelo para realizar predicciones del estado del sistema en instantes futuros es algo que se hereda del linear quadratic regulator o LQR. Esta técnica recibe su nombre de utilizar sistemas lineales con funciones de coste cuadráticas. Esta técnica propone un horizonte de eventos infinito [17]. Extender el horizonte es útil por razones de optimalidad y factibilidad, y en ausencia de restricciones, permite encontrar una solución analítica del problema. No obstante, las restricciones juegan un papel fundamental [2] en las 7 aplicaciones prácticas del control y muchos sistemas las requieren o se benefician de ellas. Es aquí donde entra en juego MPC, aplicando un horizonte finito. Esta técnica hace resoluble el problema aún con la inclusión de restricciones, utilizando métodos numéricos en vez de una solución analítica. Debe entonces encontrarse un equilibrio entre tomar 𝑁 lo suficientemente altas para acercarse lo máximo posible a la solución de horizonte infinito y lo suficientemente bajas para mantener un coste computacional adecuado. Ilustración 1. Diagrama de sistema de control genérico basado en MPC. Figura 1.2 en [1] Entre las ventajas que proporciona MPC frente a otras técnicas de control clásicas destacan las siguientes [1]: • La implementación en sistemas multivariables es una extensión del caso monovariable, con lo que no resulta especialmente compleja, a diferencia de otras técnicas de control, algunas de las cuales ni siquiera aceptan una descripción multivariable. • La adicción de restricciones es conceptualmente sencilla y puede incluirse fácilmente en el proceso de diseño. Para sistemas lineales, incluir restricciones implica pasar de una solución analítica a una numérica, generalmente de mayor costo computacional. En el caso del LQR, con horizonte infinito, la adicción de restricciones hace irresoluble al sistema. • Útil en sistemas en los que las referencias futuras son conocidas, como por ejemplo en el campo de la robótica. • Es una metodología abierta basada en unas características comunes predefinidas, lo que permite extender su aplicación a multitud de campos. Marco teórico 8 2.3 Método de la alternancia de direcciones de los multiplicadores de Lagrange El método de la alternancia de direcciones de los multiplicadores de Lagrange (ADMM) es un método para la resolución de problemas de optimización convexa. Este método surge a partir de dos métodos previos conocidos como dual ascent y método de los multiplicadores aumentados de Lagrange [16]. Así es que combina algunas de las características que hacen útiles a estos dos métodos. Por la parte de dual ascent, se tiene la descomponibilidad y por la parte de los multiplicadores, las buenas propiedades para la convergencia. Para aplicar el algoritmo se asume un problema de la forma min 𝑧,𝑣 𝑓(𝑧)+𝑔(𝑣) (5) sujeto a 𝐴𝑧+𝐵𝑣=𝑑 (6) donde 𝑓(𝑧) y 𝑔(𝑣) representan funciones convexas, 𝑧,𝑣∈ℝ𝑛 son las variables de decisión y 𝐴,𝐵∈ℝ𝑙𝑥𝑛,𝑑∈ ℝ𝑙 definen una restricción lineal para las variables de decisión [16]. Sobre este problema se puede definir un Lagrangiano aumentado 𝐿𝜌(𝑧,𝑣,𝜆)=𝑓(𝑧)+𝑔(𝑣)+𝜆𝑇(𝐴𝑧+𝐵𝑣−𝑑)+𝜌2‖𝐴𝑧+𝐵𝑣−𝑑‖22 (7) al igual que en el método de los multiplicadores. El Lagrangiano se compone de tres sumandos: la función original, las restricciones de igualdad multiplicadas por los multiplicadores de Lagrange o variable dual 𝜆 [18] y el término aumentado. Este último término consiste en la norma de la restricción multiplicada por un parámetro de penalización 𝜌 [16]. A partir de aquí, se define el algoritmo como los siguientes pasos [16]. • Condición inicial: (𝑣0,𝜆0)=(0,0) (8) 1) 𝑧𝑘+1=argmin 𝑧𝐿𝜌(𝑧,𝑣𝑘,𝜆𝑘) (9) 2) 𝑣𝑘+1=argmin 𝑣𝐿𝜌(𝑧𝑘+1,𝑣,𝜆𝑘) (10) 3) 𝜆𝑘+1=𝜆𝑘+𝜌(𝐴𝑧𝑘+1+𝐵𝑣𝑘+1−𝑑) (11) • Condición de salida: ‖𝐴𝑧𝑘+1+𝐵𝑣𝑘+1−𝑑‖2≤𝜀𝑝 (12) ‖𝑧𝑘+1−𝑧𝑘‖2≤𝜀𝑑 (13) Estos pasos siguen la misma estructura que dual ascent y el método de los multiplicadores aumentados, con un paso (o dos) de minimización de la variable (o las variables) y un paso de actualización de la variable dual. La diferencia entre ADMM y los dos anteriores es que, en caso de aplicar estos dos a un problema de la forma (5)- (6) la minimización de ambas variables 𝑧 y 𝑣 se realiza de manera simultánea, mientras que en ADMM ese paso se desdobla en dos pasos diferentes. Se puede observar en la práctica que las propiedades de convergencia de ADMM no lo hacen un algoritmo especialmente rápido, pero sí lo suficiente como para alcanzar resultados de una precisión considerable en unas decenas de iteraciones para la mayoría de aplicaciones prácticas [16]. 9 2.4 Álgebra sparsa El álgebra sparsa es un tipo de notación que permite reducir el espacio de almacenamiento y los tiempos de cómputo en cálculos relacionados con matrices de gran tamaño y gran cantidad de elementos nulos. Se trata de una herramienta matemática que está estrechamente relacionada y aporta utilidad en el cálculo numérico y en las ciencias de la computación. Una de las premisas clave para este proyecto es la eficiencia y en ese sentido es donde encuentra un lugar al uso de este tipo de notación, diseñada específicamente con la eficiencia como objetivo. Existen diversos métodos para describir matrices de forma sparsa. En general, todos consisten en guardar las componentes no nulas y los índices que las sitúan dentro de la matriz en diversos vectores. A lo largo de este proyecto se utilizan dos formatos: compressed sparse row (CSR) y compressed sparse column (CSC) [19]. Las sucesivas explicaciones se realizan sobre el primer formato, dado que el segundo es idéntico, cambiando la función realizada por las filas a las columnas. En primer lugar, debe definirse como funciona el formato y, seguidamente, se explican los dos tipos de operaciones que se han realizado a lo largo del proyecto con este tipo de matrices: el producto matriz-vector y la resolución de sistemas de ecuaciones de la forma 𝐴𝑥=𝑏. Nótese que es un tema estrechamente relacionado con la computación, así que se debe mencionar que las explicaciones dadas se basan en indexación con comienzo en el uno. Para el formato CSR, toda matriz queda descrita con un conjunto de tres vectores. El primero de ellos, vector de valores, contiene todos los valores distintos de cero, de izquierda a derecha y de arriba abajo. El segundo de ellos, vector de columnas, contiene los índices de la columna de cada uno de esos valores. Es decir, el primer elemento de este vector tiene almacenado la columna en la que se encuentra el primer elemento del vector de valores, y así sucesivamente. El último vector es el que rompe esta dinámica, el vector de filas. En este vector, la n-ésima componente del vector guarda información de la n-ésima fila. La información que guarda es en qué componente del vector de valores comienza dicha fila. El vector de filas tiene un tamaño igual al número de filas más uno. El último valor que queda almacenado en él es el número total de elementos no nulos más uno (o visto de otra forma, la longitud del vector de valores más uno). Esto cobra sentido cuando se reconstruye la matriz. El ejemplo (14) ilustra cómo funciona el formato. Los subíndices marcan la posición que ocupa el elemento en el vector o matriz correspondiente. 𝑀=(21,𝟏31,𝟐0 042,𝟐0 0 0 63,𝟑 0 0 0 0 0 0 072,𝟓0 73,𝟒83,𝟓0 0 0 104,𝟔) (14) 𝑣𝑒𝑐𝑡𝑜𝑟𝑣𝑎𝑙𝑜𝑟𝑒𝑠=(2𝟏324𝟑746𝟓768710𝟖) 𝑣𝑒𝑐𝑡𝑜𝑟𝑐𝑜𝑙𝑢𝑚𝑛𝑎𝑠=(1 2 2 5 3 4 5 6) 𝑣𝑒𝑐𝑡𝑜𝑟𝑓𝑖𝑙𝑎𝑠=(1 3 5 8 9) El formato CSR permite leer fácilmente los elementos de la matriz por filas. Para reconstruir una matriz a partir de su descripción sparsa, es necesario extraer los datos fila a fila y luego concatenar esas filas en una única matriz. Por eso, este formato resulta especialmente útil para la multiplicación matriz-vector. En ese tipo de operación, la componente n-ésima del vector resultado viene dada por el producto escalar de la n-ésima fila de la matriz por el vector multiplicado. Es así que conviene poder acceder a la matriz por filas. Para la multiplicación vector-matriz, resulta más favorable el formato CSC. Para reconstruir las filas se sigue el siguiente proceso. Se selecciona la componente del vector de filas correspondiente a la fila que se desea reconstruir como comienzo y la siguiente componente como fin. Del vector de valores, se toman las componentes de índice comprendido entre el comienzo y el fin. Así se obtienen todas las componentes no nulas de esa fila. Se repite el proceso en el vector de columnas y así se obtiene la posición exacta que ocupa cada valor dentro de la fila. En el ejemplo (15) se reconstruye la fila 2 (𝐹𝑀2) de (14). Si se quiere construir la matriz, el resto de elementos se rellenan con ceros. Si se desea hacer el producto matriz-vector, cada elemento no nulo se multiplica por el elemento correspondiente del vector, según indique el índice de columnas. Solver para MPC basado en ADMM 16 son las que definen esta condición de desigualdad, representan un lugar geométrico más sencillo: un cuadrado, cubo o equivalente en dimensiones mayores. Por tanto, las condiciones de desigualdad pueden mantenerse tal y como quedan planteadas en el MPC, resultado así el problema QP en min 𝑧12𝑧𝑇𝐻𝑧+𝑞𝑇𝑧 (34) sujeto a 𝐴𝑒𝑞𝑧=𝑏 (35) 𝑧≤𝑧≤𝑧 (36) Quedan así tres elementos entre los que establecer la conexión: las variables del problema, la función objetivo y la condición de igualdad. La manera más cómoda de hacerlo será identificando términos. Para ello, debe realizarse un paso previo, desarrollar la función (22). En primer lugar, la norma ponderada puede ser desarrollada como ‖𝑥𝑗−𝑥𝑟‖𝑄 2=(𝑥𝑗−𝑥𝑟)𝑇𝑄(𝑥𝑗−𝑥𝑟)=(𝑥𝑗𝑇𝑄−𝑥𝑟𝑇𝑄)(𝑥𝑗−𝑥𝑟) =𝑥𝑗𝑇𝑄𝑥𝑗−𝑥𝑗𝑇𝑄𝑥𝑟−𝑥𝑟𝑇𝑄𝑥𝑗+𝑥𝑟𝑇𝑄𝑥𝑟 =𝑥𝑗𝑇𝑄𝑥𝑗−2·𝑥𝑟𝑇𝑄𝑥𝑗+𝑥𝑟𝑇𝑄𝑥𝑟 (37) Debe notarse que, de los tres términos resultantes, uno de ellos solo depende de la referencia y de la matriz de ponderación. En consecuencia, este término afecta al valor final del coste, pero no a que valor de las variables de decisión lo hace mínimo y, por tanto, podemos despreciarlo del desarrollo, dejando el resultado en dos sumandos. Este mismo procedimiento puede aplicarse a las otras dos normas ponderadas que aparecen en (22). La función de coste queda entonces como 𝐽=min 𝑥,𝑢[(𝑥𝑁 𝑇𝑃𝑥𝑁−2𝑥𝑟𝑇𝑃𝑥𝑁)+∑(𝑥𝑗𝑇𝑄𝑥𝑗−2𝑥𝑟𝑇𝑄𝑥𝑗)+(𝑢𝑗𝑇𝑅𝑢𝑗−2𝑢𝑟𝑇𝑅𝑢𝑗) 𝑁−1 𝑗=0 ] (38) Si se deshacen los paréntesis, pueden reagruparse los términos de la siguiente manera. Por un lado, se agrupan los términos en los que aparece una matriz ponderada multiplicada a su derecha por un vector de variables y a su izquierda por la traspuesta de ese mismo vector. Por otro lado, se agrupan los términos que tienen una variable multiplicando a la derecha y a un vector de referencias a la izquierda. Es más, el dos que multiplica a todos los elementos de este segundo grupo puede dejarse fuera del paréntesis como factor común. El resultado es 𝐽=min 𝑥,𝑢 [(𝑥𝑁 𝑇𝑃𝑥𝑁+∑𝑥𝑗𝑇𝑄𝑥𝑗+𝑢𝑗𝑇𝑅𝑢𝑗 𝑁−1 𝑗=0 ) −2(𝑥𝑟𝑇𝑃𝑥𝑁+∑𝑥𝑟𝑇𝑄𝑥𝑗+2𝑢𝑟𝑇𝑅𝑢𝑗 𝑁−1 𝑗=0 )] (39) La forma de esta expresión es parecida a la forma de la función objetivo del problema QP, sin embargo, es necesario reagrupar más los términos, de modo que se encuentren unidas todas las variables cuadradas por un lado y las lineales por otro. Para ello, debe desarrollarse el sumatorio de la siguiente manera 𝐽=min 𝑥,𝑢[𝑥0𝑇𝑄𝑥0+⋯+𝑥𝑁−1 𝑇𝑄𝑥𝑁−1+𝑥𝑁 𝑇𝑃𝑥𝑁+𝑢0𝑇𝑅𝑢0+⋯+𝑢𝑁−1 𝑇𝑅𝑢𝑁−1 −2(𝑥𝑟𝑇𝑄𝑥0+⋯+𝑥𝑟𝑇𝑄𝑥𝑁−1+𝑥𝑟𝑇𝑃𝑥𝑁+𝑢𝑟𝑇𝑅𝑢0+⋯ +𝑢𝑟𝑇𝑅𝑢𝑁−1)] (40) 17 3.1.3 Identificación de términos En este punto, ya puede redefinirse la variable del problema, de modo que se ajuste a la del problema QP. Simplemente se construye un nuevo vector que contenga a todas las variables, los estados y señales de control para cada instante. No existe una forma única de distribuir las variables dentro del vector. Para este proyecto, se ha seleccionado la siguiente 𝑧= ( 𝑢0 𝑥1 𝑢1 𝑥2 𝑢2 ⋮ 𝑥𝑁−1 𝑢𝑁−1 𝑥𝑁 ) (41) Esta forma hace más fácil la construcción del resto de elementos. Además, las variables que interesan desde el punto de vista del control (las señales de control en el instante actual) ocupan las primeras posiciones del vector. Si se reescriben todos los sumandos en (40) como una operación matricial, de modo que las variables de la ecuación resulten en (41), las matrices de ponderación formarán una nueva matriz. Esta matriz es el equivalente de la matriz 𝐻 del problema QP. Por tanto, la formación de esta matriz queda como 𝐻= ( 𝑅𝑄𝑅𝑄𝑅⋱𝑄𝑅𝑃 ) (42) Para los términos lineales, puede aplicarse el mismo razonamiento, siendo el resultado 𝑞=− ( 𝑅𝑢𝑟 𝑄𝑥𝑟 𝑅𝑢𝑟 ⋮ 𝑄𝑥𝑟 𝑅𝑢𝑟 𝑃𝑥𝑟 ) (43) Cabe comentar un par de detalles. En primer lugar, mientras el término lineal del MPC aparecía restando, en la formulación QP aparece sumando, por lo que todas las componentes del vector aparecen con el signo cambiado. En segundo lugar, al ser las matrices ponderadas simétricas por definición, es equivalente multiplicarles el vector traspuesto de referencia por la izquierda, que multiplicarles el mismo por la derecha sin trasponer. Además, toda la expresión puede dividirse entre dos sin alterar la solución óptima. Esto se debe a que la solución óptima que se busca es la variable de decisión óptima, y no el valor óptimo de la función (que sí quedaría alterado). Una vez se ha formado el vector de variables, atendiendo a (36), los vectores de cota superior e inferior resultan en Solver para MPC basado en ADMM 18 𝑧= ( 𝑢𝑥⋮𝑢𝑥 ) 𝑧= ( 𝑢𝑥⋮𝑢𝑥 ) (44) El único elemento restante para terminar de adaptar MPC a QP son las restricciones de igualdad. En el problema planteado, las restricciones de igualdad vienen dadas por la condición inicial del problema (25). Si se sustituye esta expresión en el modelo de predicción se obtiene 𝑥1=𝐴𝑥0+𝐵𝑢0⁡→⁡𝑥1=𝐴𝑥(𝑘)+𝐵𝑢0 (45) Redistribuyendo los miembros de la ecuación, con los datos numéricos en el primer miembro y las variables en el segundo −𝐴𝑥(𝑘)=𝐵𝑢0−𝑥1 (46) Esta ecuación puede construirse para todos los instantes del modelo de predicción hasta el horizonte de eventos. Para el segundo instante, la ecuación es 𝑥2=𝐴𝑥1+𝐵𝑢1⁡→⁡𝐴𝑥1+𝐵𝑢1−𝑥2=0 (47) Este proceso puede aplicarse a todos los instantes, hasta 𝐴𝑥𝑁−1+𝐵𝑢𝑁−1−𝑥𝑁=0 (48) de modo que el conjunto resultante de ecuaciones sigue un patrón. Ese sistema de ecuaciones puede escribirse de forma matricial como 𝐴𝑒𝑞𝑧=𝑏. Es exactamente la estructura de las restricciones de igualdad. Por tanto, la matriz y el vector independiente se construyen como 𝐴𝑒𝑞= ( 𝐵 −𝐼 0 0 𝐴 𝐵 0000 0 0 −𝐼 0 0 𝐴 𝐵 −𝐼 ⋱000 −𝐼 0 0 𝐴 𝐵 −𝐼 ) (49) 𝑏= ( −𝐴𝑥(𝑘) 00⋮0 ) (50) Con ello, ya se encuentran redefinidos todos los objetos matemáticos que conforman el problema QP en base a los valores, variables, parámetros y datos del MPC. El solver a elaborar será el encargado de resolver el problema QP y estas matrices y vectores serán los datos o ingredientes que deben proporcionárseles para definir el problema concreto que se esta resolviendo. Pues una vez definidos los ingredientes, es momento de ver las instrucciones que el solver debe ejecutar para encontrar la solución óptima del problema. 19 3.2 Uso de ADMM para resolver MPC Hasta aquí se han desarrollado las matemáticas necesarias para describir la técnica de control MPC con la formulación típica de los problemas de optimización QP. Ahora, es necesario resolver ese problema QP y para ello se ha contado con el algoritmo ADMM. Este es un algoritmo genérico para resolución de distintos tipos de problemas de optimización. Como su alcance abarca más allá de problemas QP, la función objetivo general suele venir descrita como (5)-(6). Esta función (y sus restricciones) no se corresponde con la formulación usada para MPC una vez adaptada a QP. Por tanto, es necesario hacer algunos ajustes más a la descripción del problema a la que se ha llegado en el apartado anterior antes de aplicar ADMM. Así mismo, las instrucciones generales que conforman el algoritmo (8)-(13) deben adaptarse a la formulación que está siendo usada. 3.2.1 Adaptación de la función objetivo En primer lugar, resalta que la función objetivo propuesta (5) para ADMM está compuesta a su vez por dos funciones distintas, cada una de ellas dependiente de una variable distinta. En consecuencia, la variable del problema QP debe romperse en dos variables distintas y los elementos que componen la función objetivo deben reagruparse en dos funciones individuales. Para la restricción propuesta en el enunciado del problema genérico para ADMM (6) se seleccionan los valores de las matrices y el término independiente de forma que se pueda reescribir la ecuación como min 𝑧,𝑣 12𝑧𝑇𝐻𝑧+𝑞𝑇𝑧+𝔗𝑒𝑞(𝑧)+𝔗𝑖𝑛𝑒𝑞(𝑣) (51) sujeto a 𝑧=𝑣 (52) Así queda el problema original en dos subproblemas desacoplados a excepción de una ecuación lineal (52) que los relaciona. Sin embargo, deja la incógnita de como expresar las restricciones de la formulación QP (35)-(36). Estas restricciones se agregan como funciones indicadoras 𝔗𝑒𝑞(𝑧) y 𝔗𝑖𝑛𝑒𝑞(𝑣) a la función objetivo, del modo que se expresa a continuación. La función objetivo queda entonces de la siguiente manera (51). Esta nueva función objetivo está compuesta por dos funciones individuales en variables distintas (requerimiento de la formulación para ADMM considerada). Las funciones indicadoras siguen una determinada estructura. Se trata de funciones a trozos que solo toman dos valores: cero e infinito. Valen cero cuando se cumplen la restricción que las generan (es decir, cuando la variable pertenece al conjunto que define la ecuación de la restricción) e infinito en caso contrario. De este modo, mientras que se cumplan las restricciones, las funciones indicadoras, que son un sumando más de la función objetivo, no aportan nada. Esto es lo deseado, ya que mientras que se cumplan las restricciones, la función objetivo se mantiene sin alterar. En caso de que no se cumpla la restricción, se está añadiendo un sumando de valor infinito a la función objetivo, de modo que ese valor de las variables que incumple las restricciones jamás será el valor que hace mínima la función. Las funciones indicadoras introducidas en el problema tienen la estructura 𝔗𝑒𝑞={0 +∞ si⁡𝐴𝑒𝑞𝑧=𝑏 en⁡otro⁡caso (53) 𝔗𝑖𝑛𝑒𝑞={0 +∞ si⁡𝑧≤𝑣≤𝑧 en⁡otro⁡caso (54) Volviendo a (51), la función de la variable z (tener en cuenta que ambas variables son iguales, es decir, la variable de decisión), contiene tres términos. El primer de ellos es el término cuadrático de la formulación QP. El segundo de ellos es el término lineal de la formulación original. El tercer término es la función indicadora que representa Solver para MPC basado en ADMM 20 a la restricción de igualdad. Por otro lado, la función de la variable v, incluye solo a la función indicadora que representa las restricciones de desigualdad. En el apartado 2.3 se comentó como el método ADMM requiere el uso del lagraniano. El Lagrangiano resultante de esta nueva reformulación del problema es 𝐿𝜌(𝑧,𝑣,𝜆)=12𝑧𝑇𝐻𝑧+𝑞𝑇𝑧+𝔗𝑒𝑞(𝑧)+𝔗𝑖𝑛𝑒𝑞(𝑣)+𝜌2‖𝑧−𝑣‖2+〈𝜆,𝑧−𝑣〉 (55) 3.2.2 Desarrollo del algoritmo Una vez el problema ha sido adaptado a la expresión deseada, es momento de aplicar el algoritmo ADMM. De forma general, sigue la siguiente estructura (8)-(13). A continuación, se analiza cada una de las partes que lo componen de forma individual. La condición inicial se mantiene igual, inicializando a cero los vectores correspondientes. Lo siguiente son los tres pasos del algoritmo. El primero corresponde a la variable z y el segundo a la variable v. En cada uno de ellos, se halla el valor de la variable que minimiza el Lagrangiano sujeto a los valores de esa iteración para las otras. Es decir, la función a minimizar es dependiente solo de la variable que se pretende calcular. Debido a eso, en cada paso se toma como función objetivo solo a la función individual para esa variable más el par de sumandos que se añaden con el Lagrangiano, dando sentido a por qué separar la función objetivo en dos funciones diferentes. El resultado entonces para el primer paso del algoritmo es 1) 𝑧𝑘+1=argmin 𝑧12𝑧𝑇𝐻𝑧+𝑞𝑇𝑧+𝜌2‖𝑧−𝑣𝑘‖2+〈(𝜆𝑘)𝑇,𝑧〉 (56) sujeto a 𝐴𝑒𝑞𝑧=𝑏 Puede observarse que la expresión vuelve a ser un problema de optimización en una única variable vectorial y, por tanto, puede deshacerse la función indicadora para ser expresada como una restricción de igualdad de nuevo. Si se desarrolla la norma y se agrupan términos, puede llegarse a la expresión 12𝑧𝑇𝐻𝑧+𝑞𝑇𝑧+(𝜆𝑘)𝑇𝑧+𝜌2𝑧𝑇𝐼𝑧−𝜌(𝑣𝑘)𝑇𝑧+𝜌2(𝑣𝑘)𝑇𝑣𝑘 =12𝑧𝑇(𝐻+𝜌𝐼)𝑧+(𝑞+𝜆𝑘−𝜌𝑣𝑘)𝑇𝑧=12𝑧𝑇𝐻𝑧𝑧+(𝑞𝑧𝑘)𝑇𝑧 (57) Queda entonces el siguiente problema 𝑧𝑘+1=argmin 𝑧12𝑧𝑇𝐻𝑧𝑧+〈(𝑞𝑧𝑘)𝑇,𝑧〉 con⁡𝑞𝑧𝑘=(𝑞+𝜆𝑘−𝜌𝑣𝑘) (58) sujeto a 𝐴𝑒𝑞𝑧=𝑏 Intuitivamente se observa una estructura parecida con el problema QP originalmente planteado. La diferencia, más allá del valor de los datos, es la ausencia de restricciones de desigualdad. Resulta que en esta situación sí existe una solución analítica explicita al problema [18] dada por 𝑊𝜇=−(𝐴𝑒𝑞𝐻𝑧−1𝑞𝑧𝑘+𝑏) con⁡𝑊=(𝐴𝑒𝑞𝐻𝑧−1𝐴𝑒𝑞 𝑇) (59) 𝑧𝑘+1=−𝐻𝑧−1(𝐴𝑒𝑞 𝑇𝜇+𝑞𝑧𝑘) (60) 21 Aplicando la solución, el cálculo de 𝑧𝑘+1 queda expresado como una única expresión que puede escribirse fácilmente en código. Para el segundo paso, puede aplicarse un desarrollo idéntico al primero, con lo que 2) 𝑣𝑘+1=argmin 𝑣𝜌2‖𝑧𝑘+1−𝑣‖2−〈(𝜆𝑘)𝑇,𝑣〉 (61) sujeto a 𝑧≤𝑣≤𝑧 pasa a ser 𝑣𝑘+1=argmin 𝑧12𝑣𝑇(𝜌𝐼)𝑣+〈(𝑞𝑣𝑘)𝑇,𝑣〉 con⁡𝑞𝑣𝑘=(−𝜆𝑘−𝜌𝑧𝑘+1) (62) sujeto a 𝑧≤𝑣≤𝑧 siguiendo 𝜌2𝑣𝑇𝐼𝑣−𝜌(𝑧𝑘+1)𝑇𝑣+𝜌2(𝑧𝑘+1)𝑇𝑧𝑘+1⁡−(𝜆𝑘)𝑇𝑣 =12𝑣𝑇(𝜌𝐼)𝑣+(−𝜆𝑘−𝜌𝑧𝑘+1)𝑇𝑣=12𝑣𝑇(𝜌𝐼)𝑣+(𝑞𝑣𝑘)𝑇𝑣 (63) En este caso, la particularidad del problema resultante es que se encuentra totalmente desacoplado. Esto se debe a la forma de la restricción y a que la norma no introduce acoplamientos. Debido al desacople, el problema cuadrático en una variable vectorial (varias variables) puede descomponerse en un número de problemas cuadráticos de una variable escalar (una única variable) independientes entre sí. El número de problemas independientes es igual al tamaño de v. Los problemas cuadráticos en una única variable describen una parábola en el plano. Es conocimiento matemático básico que el extremo relativo de una parábola se encuentra en el vértice. Las restricciones de desigualdad representan un segmento en el que se observa la parábola, un intervalo en el que debe encontrarse la solución. Por tanto, el mínimo absoluto debe ser el vértice o uno de los límites del intervalo. El procedimiento queda entonces como sigue 𝑣𝑖𝑘+1=max[𝑧𝑖,min[𝑧𝑖,−(𝑞𝑣𝑘)𝑖 𝜌]] (64) Se comprueba que valor es menor, si el vértice o la cota superior. Después, se toma el valor mayor entre la cota inferior y el resultado de esta operación. Así, si el vértice se encuentra entre las dos cotas, el resultado es el vértice. Si el vértice se encuentra por encima de la cota superior, el resultado será la cota superior, y si se encuentra por debajo de la cota inferior, este será el resultado final. El tercer y último paso de ADMM es (11), que tras aplicar la restricción planteada en (52), resulta en 3) 𝜆𝑘+1=𝜆𝑘+𝜌(𝑧𝑘+1−𝑣𝑘+1) (65) Este paso no requiere mayor desarrollo, pues en su estado actual ya es implementable como código. Finalmente, como fin del algoritmo, deben cumplirse unas condiciones de terminación o salida. La primera, (12), asegura que el resultado óptimo encontrado cumple las restricciones propuestas. En el enunciado, formulado de manera analítica, se esperaría que la expresión que representa las restricciones fuese igual a cero como condición de salida. Sin embargo, al tratarse de un método numérico, es suficiente con que dicha expresión Solver para MPC basado en ADMM 22 tome un valor menor a una determinada tolerancia, que puede ajustarse como parámetro de diseño del solver. La segunda condición, (13), permite tener certeza de que se ha hallado un valor estable. Aplica el mismo razonamiento sobre las tolerancias. Cabe añadir una tercera condición, de carácter computacional y no relacionada con el problema matemático, que evite que el algoritmo se trabe en un bucle infinito. El resultado de aplicar estas condiciones al problema exacto planteado es ‖𝑧𝑘+1−𝑣𝑘+1‖2≤𝜀𝑝 (66) ‖𝑧𝑘+1−𝑧𝑘‖2≤𝜀𝑑 (67) 𝑘≥𝑘𝑚𝑎𝑥 (68) Una vez que se cumplen esas condiciones, el algoritmo ha finalizado. El valor de la variable 𝑧𝑘+1 es el valor óptimo para el que la función objetivo es mínima. Por la construcción del vector de variables mostrada en el apartado anterior, las primeras 𝑚 componentes, siendo 𝑚 el número de señales de control, son la acción de control que se debe aplicar al sistema durante ese tiempo de muestreo para que se alcance la referencia. 3.3 Álgebra sparsa aplicada a ADMM El desarrollo de solvers para problemas de optimización es un campo sobradamente estudiado. Como aporte de este proyecto, se ha pretendido tomar esos conocimientos y aplicarlos al desarrollo en un lenguaje novedoso como es Julia. En ese sentido se tiene que la eficiencia es un aspecto central de este proyecto, tal y como se ha comentado en los apartados introductorios, y por ello no se ha limitado su búsqueda a cuestiones de lenguaje. Desde el punto de vista del propio desarrollo, también pueden tomarse medidas para reducir el coste computacional de todo el aparato matemático que se está manejando. Durante la aplicación del algoritmo ADMM, se realizan numerosas operaciones con vectores y matrices. En base al número de estados y de señales de control del sistema al que se enfrente el solver, estas entidades algebraicas pueden llegar a tener unas dimensiones considerables. Es más, se da el caso de que algunas de estas matrices cuenten con gran cantidad de componentes nulas. Así es el caso de la matriz 𝐻, que es diagonal a bloques. Estos espacios nulos ocupan gran cantidad de memoria y suponen una serie de cálculos que no proporciona información útil sobre el problema. Para ahorrar esfuerzo computacional en la ejecución del algoritmo ADMM se ha recurrido a la notación sparsa en puntos claves del mismo. El principal conjunto de operaciones matriciales que realizan en el desarrollo de ADMM es (59)-(60). El par de ecuaciones puede ser expresado como 𝑊𝜇=𝑏 (69) con⁡𝑏=𝐴𝑒𝑞𝐻𝑧−1𝑞𝑧𝑘+𝑏 −𝐻𝑧𝑧𝑘+1=𝜇 (70) con⁡𝜇=𝐴𝑒𝑞 𝑇𝜇+𝑞𝑧𝑘 Puede observarse que, en realidad, se trata de un par de sistemas de ecuaciones de la forma 𝐴𝑥=𝑏. Además, los términos independientes están formados por una serie de operaciones que pueden reducirse a un producto matriz-vector al que se le suma otro vector. La inversa de 𝐻𝑧, al ser esta diagonal a bloques y definida positiva, sigue la misma estructura, por lo que convertirla a CSR es muy útil para resolver (70). La descomposición LDL de 𝑊 también es muy sparsa y útil para resolver (69). Por ello, tanto la pareja de productos matriz-vector como ambos sistemas de ecuaciones pueden resolverse de manera sparsa. Esto conlleva que los datos que se proporcionen al solver, tanto como las operaciones realizadas, requieran modificarse de cara a implementar este par de ecuaciones. 23 3.4 Implementación de un controlador MPC en bucle cerrado De todo este conjunto de operaciones es de lo que se componen las entrañas del controlador. Pero hablar de un controlador aislado tiene poco sentido. Es necesario colocarlo en su contexto y ver como se relaciona con los demás participantes en el sistema. En especial, uno de los conceptos más importantes a tener en cuenta en la estructura del sistema es la realimentación, independientemente de la técnica de control utilizada. La realimentación es omnipresente en los sistemas de control, pues permite al sistema controlado adaptarse a cambios en la referencia o rechazar los efectos de las perturbaciones. Para este proyecto, al estar centrado en la técnica de control principalmente, se ha decidido no implantar una estructura de control muy complicada. Como ya se ha mencionado, se ha descartado observar las salidas del sistema y lo que se pretende controlar son los estados del sistema. Se ha establecido que el sistema de control está compuesto simplemente por el controlador que implementa MPC y por el sistema físico a controlar. El controlador recibe como entradas la referencia (estados y señales de control) y calcula la señal de control adecuada, la óptima, obtenida como solución del problema QP. El controlador entrega las señales de control al sistema, que modifica su estado en consecuencia. La única interacción más que aparece en el sistema es la realimentación. Esta es representada por un flujo de información del sistema al controlador. El controlador por MPC debe ser capaz de leer el estado actual del sistema en cada tiempo de muestreo debido a que lo necesita para realizar las operaciones del cálculo de la señal de control óptima. Ilustración 2: Diagrama de control en bucle cerrado para MPC Para poder testear el solver en funcionamiento, se debe realizar una simulación del sistema. Para ello, se enfrenta al controlador a un modelo del sistema. Partiendo de unos valores iniciales, se simulan a través de un bucle los ciclos de muestreo. En cada tiempo de muestro, actualizan los valores de 𝑥(𝑘) y (𝑥𝑟,𝑢𝑟) en primer lugar. Luego se resuelve el MPC. A continuación, se simula el sistema realizando las operaciones del modelo. No es estrictamente necesario para el funcionamiento del sistema, pero, finalmente, puede guardarse los valores de las variables que se consideren, para su posterior análisis. La Tabla 1 enseña el pseudocódigo que ejemplifica este funcionamiento. Tabla 1. Pseudocódigo de simulación 1 𝑥(𝑘)=𝑥0 2 𝐟𝐨𝐫⁡1⁡𝐭𝐨⁡n_iteraciones 3 Actualizar⁡(𝑥𝑟,𝑢𝑟) 4 Resolver⁡MPC⁡→𝑢0∗ 5 Simular⁡sistema⁡𝑥(𝑘+1)=𝐴𝑥(𝑘)+𝐵𝑢0∗ 6 (Guardar⁡valores⁡𝑢(𝑘),𝑧∗,𝑘) 7 𝐞𝐧𝐝⁡𝐟𝐨𝐫 25 4 DESARROLLO DEL SOLVER El solver constituye el elemento central del proyecto. Aunque el objetivo final se limite a crear el solver y comprobar su correcto funcionamiento en simulación, el enfoque desde el que se realiza es que fuese ejecutable, de manera eficiente, por un dispositivo electrónico capaz de controlar un sistema físico. Se hace especial énfasis en la eficiencia para que la gama de dispositivos que pudiesen ejecutar el código abarque los dispositivos más modestos posibles. Esta característica marca el diseño del mismo desde su propia estructura, que es presentada a continuación. En este apartado se realiza una descripción del conjunto de códigos que conforman la totalidad del solver, la estructura que siguen y su diseño. Se explican las decisiones de diseño tomadas y las dificultades afrontadas durante el desarrollo. De las explicaciones se deducen y se exponen los pseudocódigos de las funciones implementadas. El código final implementado en los respectivos lenguajes se adjunta en los anexos de este documento. 4.1 Desarrollo en Matlab Una de las decisiones de diseño más importantes es el lenguaje en el que debe ser implementado. Como ya se ha explicado, la elección tomada fue Julia. Sin embargo, debido a la cantidad de herramientas que ofrece y la potencia de su entorno de desarrollo, se decidió hacer una versión preliminar en Matlab. Esto ha permitido familiarizarse de forma más cómoda con la idiosincrasia de MPC. Una vez escrito todo el código en Matlab, ha sido transcrito a Julia haciendo los ajustes correspondientes. 4.1.1 Descripción del solver Un solver es un programa que aplica una serie de operaciones matemáticas sobre un conjunto de datos a través de un algoritmo para proporcionar la solución a un determinado problema matemático caracterizado por ese conjunto datos. En el contexto de este proyecto ese problema matemático es un problema de optimización tipo QP. Como el problema de optimización concreto deriva de la descripción del sistema de control MPC, recontextualizado, implica que es necesario distinguir entre dos bloques diferenciados entre sí para la implementación: la adecuación del problema de control al problema matemático y la solución del problema matemático en sí. Esta es la principal diferenciación que marca la estructura del código desarrollado. Esto hace que el código se presente como dos funciones principales. Una de ellas es el solver en sí y la otra, la función que genera los datos o ingredientes del solver a partir de los datos del sistema. Esta idea de dividir el problema en dos tiene más calado si cabe a la hora de plantear una situación real. Para un controlador real que implemente MPC, su trabajo para cada tiempo de muestreo es resolver el problema QP, es decir, ejecutar el solver. Es así como se calcula la señal de control que debe implementarse en el siguiente ciclo de muestreo. Todo el tema relacionado con la adaptación del modelo al solver es algo relacionado con la caracterización del sistema, previo a la puesta en marcha del controlador. Por tanto, resulta ineficiente realizar esos cálculos siempre que se ejecute el solver y basta con realizarlos durante la inicialización del sistema. Es más, tratándose de controladores con pocos recursos, los ingredientes pueden calcularse de modo externo al controlador y quedar almacenados en una memoria interna, de modo que el solver simplemente pueda leerlos desde ahí cada vez que se ejecute. A continuación, se describen cada una de las partes que componen el solver. El grueso de las funciones consiste en ejecutar en forma de instrucciones de código los cálculos de las ecuaciones presentadas en el Capítulo 3. Sin embargo, debido a las particularidades de los lenguajes y la computación, conviene realizar pasos intermedios, preparar los datos y arreglar las operaciones a la forma más conveniente para cada lenguaje. Por tanto, en los Desarrollo del solver 32 4.2 Desarrollo en Julia Desarrollar el solver suponía enfrentarse al problema MPC, familiarizarse con él y entenderlo en profundidad. A raíz de eso, se eligió Matlab como entorno de desarrollo, pues suponía una potente herramienta de cara al aprendizaje, más consolidada en cuestiones de control y simulación de sistemas frente al joven y novedoso Julia. Son mayor la cantidad de recursos que facilita para trabajar con cálculo numérico y sistemas de control, además de tener un ecosistema más autocontenido. A pesar de ello, este proyecto está orientado a la implementación en Julia, que es lo que suponía un reto innovador. Tabla 8. Pseudocódigo del algoritmo QDLDL 𝐄𝐧𝐭𝐫𝐚𝐝𝐚:⁡𝐷𝑖𝑛𝑣,𝐿𝐶𝑆𝐶⁡(sparsa∗),𝑏 1 𝑠𝑜𝑙=𝑏 2 𝐟𝐨𝐫⁡i=1⁡𝐭𝐨⁡dimension_D 3 inicio=vector_columna_Li 4 fin=vector_fila_Li+1 5 𝐟𝐨𝐫⁡j=inicio⁡𝐭𝐨⁡fin 6 𝑠𝑜𝑙vector_filas_Lj=𝑠𝑜𝑙vector_filas_Lj+vector_valores_Lj·𝑠𝑜𝑙𝑖 7 𝐞𝐧𝐝⁡𝐟𝐨𝐫 8 𝐞𝐧𝐝⁡𝐟𝐨𝐫 9 𝐟𝐨𝐫⁡i=1⁡𝐭𝐨⁡dimension_D 10 𝑠𝑜𝑙𝑖=𝑠𝑜𝑙𝑖·𝐷𝑖𝑖𝑛𝑣 11 𝐞𝐧𝐝⁡𝐟𝐨𝐫 12 𝐞𝐧𝐝⁡𝐟𝐨𝐫 13 𝐟𝐨𝐫⁡i=1⁡𝐭𝐨⁡dimension_D 14 inicio=vector_columna_Li 15 fin=vector_fila_Li+1 16 𝐟𝐨𝐫⁡j=inicio⁡𝐭𝐨⁡fin 17 𝑠𝑜𝑙𝑖=𝑠𝑜𝑙𝑖+vector_valores_Lj·𝑠𝑜𝑙vector_filas_Lj 18 𝐞𝐧𝐝⁡𝐟𝐨𝐫 19 𝐞𝐧𝐝⁡𝐟𝐨𝐫 𝐒𝐚𝐥𝐢𝐝𝐚:⁡𝑠𝑜𝑙 ∗Se⁡supone⁡que⁡𝐿𝐶𝑆𝐶⁡es⁡una⁡matriz⁡en⁡notación⁡CSC⁡ compuesta⁡por:vector_valores,vector_columnas,vector_filas Una vez que el solver ha sido desarrollado conceptualmente, implementarlo en un lenguaje u otro poco tiene que ver con el solver en sí y más con el propio lenguaje. Por tanto, desarrollar el solver en Julia consiste mayoritariamente en transcribir el código escrito en Matlab. Las explicaciones del funcionamiento del solver ya se han hecho con Matlab y no es necesario repetirlas. En este apartado se comentan las necesidades específicas que se han tenido que abordar para adaptar el código de un lenguaje a otro. Además, en el apartado anterior se han comentado estrictamente las distintas partes que componen el solver y cómo funcionan cada una de ellas. En este apartado se realiza una explicación más amplia de la estructura en la que se encuentra trabajando el solver. Esto se debe a que, aunque el solver se haya transcrito e implementado en Julia, parte del trabajo auxiliar que se requiere sigue siendo ejecutado en Matlab. Esto está relacionado con el hecho varias veces reiterado de que la función generatriz no tiene por qué ser ejecutada por el dispositivo que ejecuta el solver. 33 4.2.1 Estructura de la solución propuesta en Julia La idea es que todos los ingredientes y los parámetros respectivos al sistema sean generados desde un script de Matlab. El script crea un fichero con todos estos datos, que deben ser leídos por el solver. Este fichero podría ser cargado en la memoria del dispositivo que ejecutase el solver. En un entorno real, la función que ejecuta el algoritmo del solver no tiene capacidad para leer los datos directamente de la memoria. Se necesitaría un script de adaptación (o bien adaptar la función solver para dotarla de esa capacidad) que se encargase de leer los datos del fichero y, en cada tiempo de muestreo, recibir los valores actuales de estado y referencia, para entregárselo todo al solver. El solver, al ser una función llamada desde ese fichero, le devuelve el resultado de la acción de control. Debería ser el script el encargado de mandar al controlador para que este aplique acción de control al sistema. Esto es un caso hipotético, el fin último del solver. Sin embargo, el alcance de este proyecto se ha limitado a que el solver realice los cálculos correctamente, por lo que las pruebas se han realizado en un entorno simulado. En este caso, es un propio script de Julia el encargado de simular el sistema físico, y es el que realiza las tareas de leer el fichero y entregarle los datos correspondientes al solver. Ilustración 3. Estructura de la solución propuesta en Julia Por lo tanto, la función generatriz, uno de los dos principales bloques del programa, no necesita ser transcrita a Julia. Llamar a esta función, junto a otras tareas, es el trabajo del script de Matlab. Entre esas otras tareas está definir los parámetros del modelo y los de ADMM o calcular el estado de referencia, principalmente. Los datos Desarrollo del solver 34 generados con el cálculo del solver también son recogidos y almacenados en un fichero de Matlab. Estos datos incluyen las señales de control óptimas, las variables óptimas resultado de la resolución del problema QP y el número de iteraciones en calcular el resultado. Cuando el solver ha terminado su funcionamiento es cuando se genera este fichero. Para trabajar con archivos de Matlab desde Julia se ha utilizado el paquete MAT. Los resultados pueden ser visualizados posteriormente en Matlab para realizar depuración del funcionamiento del solver. Un esquema de esta estructura se representa en Ilustración 3. 4.2.2 Particularidades de Julia En la transcripción del código del solver a Julia hay que enfrentarse a algunas particularidades de los lenguajes. Para modular lo máximo posible el trabajo en Matlab se ha recurrido al uso de estructuras. Como se manejan una gran variedad de variables, conviene agruparlas según a las distintas partes del problema que hacen referencia. Se reparten las variables en grupos como parámetros del modelo, parámetros de ADMM o ingredientes. Además de mantener el código limpio y ordenado, esto es útil a la hora de trabajar con las distintas funciones creadas, pues cada una de ellas suele relacionarse con partes concretas del problema. Por ejemplo, la función generatriz tiene poca relación con los parámetros del ADMM. Por tanto, en lugar de recibir como entrada un sinfín de variables diferentes, recibe como entrada el modelo del sistema y el parámetro 𝜌 de forma independiente. En la transcripción a Julia, se ha descartado el uso de estructuras por dos motivos: el uso de MAT y que requieren ser predefinidas. Mientras que en Matlab las estructuras son definidas de manera dinámica cuando se les asignan valores, en Julia estas necesitan ser declaradas con antelación. Además, los datos almacenados en ficheros de Matlab como estructuras son interpretados por la librería MAT como diccionarios al ser leídos. Es así que se han sustituido todas las estructuras por diccionarios en Julia. Otro fenómeno de la idiosincrasia de Julia a tener en cuenta es el local scope [21]. Julia, como una gran cantidad de lenguajes modernos, cuenta con declaración implícita de variables en la asignación de valores. Esto quiere decir que, al asignar un valor a una variable, si esa variable no se encuentra previamente declarada, se declara de manera automática en esa instrucción. Pues bien, en Julia se definen dos tipos de scope o ámbito de la variable, global y local, como en la mayoría de lenguajes. Además, dentro del ámbito local, se diferencia entre soft scope y hard scope. El primero es el que prevalece en los bucles, que es el que interesa para esta cuestión, pues en cada ejecución del solver se ejecuta un bucle para el cálculo del valor óptimo. Cuando se realiza una declaración implícita dentro del bucle, a la variable declarada se le asigna un ámbito local. Con ámbito local se restringe su uso al interior del bucle en el que fue declarada. Esta restricción es llamada block scope, es decir, el ámbito local no abarca la función completa, sino únicamente el bloque (bucle en este caso) donde fue declarada la variable. Si en el mismo código se utiliza otra variable del mismo nombre fuera del bucle, el soft scope indica que lo siguiente es lo que debe suceder. Si esa variable está declarada fuera del bucle, Julia envía un mensaje de warning de ambigüedad, indicando que se está declarando una variable local dentro del bucle con el mismo nombre que una variable ajena al bucle. Esto es útil para un programador inexperto en el uso de Julia, que puede no saber que la variable que está siendo asignada en el interior del bucle es distinta a la variable externa. En caso de que la variable no esté declarada fuera del bucle, salta un error indicando que esa variable no existe. Esto es así pues una vez terminado de ejecutar el bucle, la variable local declarada implícitamente deja de existir, por lo que si se intenta asignar su valor a otra variable fuera del bucle, aparece un error. Este último ejemplo es el caso que ocupa a este proyecto. Dentro del bucle se realizan los cálculos pertinentes para obtener el resultado. La mayoría de las variables usadas son circunstanciales, requeridas para las operaciones, pero su valor final resulta indiferente. Estas no suponen ningún problema. Las variables que sí almacenan datos numéricos útiles deben ser entregadas como resultado de la función fuera del bucle. Por ello, deben ser declaradas como globales, para que esa información no desaparezca al finalizar los cálculos. 35 5 RESULTADOS NÚMERICOS Con el solver desarrollado, es necesario realizar unos experimentos prácticos para poner en valor la calidad de su funcionamiento. Para ello, se simula el controlador frente a dos sistemas distintos. El primero de ellos, compuesto por cuatro tanques de agua, y el segundo, por una serie de masas conectadas a través de muelles. En este capítulo se otorga al lector una descripción matemática detallada de los mismos. Para comparar los resultados, se realizan las simulaciones con tres solvers diferentes. Dos de ellos son el solver desarrollado en este proyecto, pero uno es simulado en Matlab mientras que el otro es simulado en Julia. El tercer solver es COSMO, un solver puntero en el estado del arte, desarrollado en la Universidad de Oxford, capaz de obtener resultados muy eficientes y que cuenta con un montón de funcionalidades. 5.1 «Conic operator splitting method for convex conic problems» El solver COSMO («COSMO: A conic operator splitting method for convex conic problems») es un solver desarrollado por el University of Oxford Control Group, basado en el método de resolución de problemas de optimización convexa del mismo nombre [22]. Este solver está implementado en el lenguaje Julia, por lo que supone una referencia ideal para medir la eficacia del solver ADMM programado. Se trata de un solver de propósito más general, enfocado en la resolución de todo tipo de problemas cuadráticos convexos y cónicos, es decir, problemas que siguen la formulación min 𝑥12𝑥𝑇𝑃𝑥+𝑞𝑇𝑥 (76) sujeto a 𝐴𝑥+𝑠=𝑏 (77) 𝑠∈𝒦 (78) con 𝒦 siendo una composición de conos y conjuntos convexos. El problema QP descrito durante este proyecto puede ser descrito en base a estas ecuaciones, por lo que COSMO puede utilizarse como solver para MPC. Entre sus ventajas figuran cualidades como versatilidad, uso de aceleradores, detección de inviabilidad, chordal decomposition o warm starting. Es importante recalcar que se trata de un solver muy potente, con muchas adicciones que permiten acelerar la resolución del problema, incluso para problemas más complejos y en mayores dimensiones del que se pretende resolver aquí. Por tanto, debe considerarse que no se pretenden mejorar los tiempos y ejecuciones de COSMO, sino usarlos como referencia del nivel de eficacia al que se encuentran los solvers en la vanguardia actual y así contar con una medición de dónde se encontraría el solver desarrollado dentro de este contexto. Al igual que en el propio solver ADMM, que funcionaba como algoritmo para resolver problemas QP, COSMO es una herramienta para optimización, en este caso incluso más general. Por ello debe ser adaptada a la idiosincrasia del problema MPC. Dentro de la estructura de Ilustración 3, COSMO ocupa la posición de solver. Esto significa que requiere del resto de la estructura señalada para adecuarse al control MPC. Además, a pesar de reutilizar la estructura diseñada, esta debe adecuarse a las particularidades de COSMO. En primer lugar, como ya se ha mencionado, es necesario ajustar el problema convexo génerico que resuelve COSMO al problema QP que se está resolviendo en este caso. En un principio, en cuanto a la función objetivo, simplemente debe señalarse el resultado obvio de que la variable de decisión en (31) es 𝑧 y la matriz Hessiana se Resultados númericos 36 denota como 𝐻 y, por tanto, 𝑥=𝑧 y 𝑃=𝐻. La auténtica diferencia la marcan las restricciones. Para adaptar a las restricciones (35)-(36) la formulación de (77)-(78) se aplica (𝐶𝐼)𝑧+(𝑠1 𝑠2)=(𝑑0) (79) (80) con⁡𝑠1={0}⁡,𝒦={𝑧∶⁡𝑢≤𝑧≤𝑙}⇒𝑢≤𝑠2≤𝑙 teniendo en cuenta que (35) se ha reescrito como 𝐶𝑧=𝑑. De (36) y (80) se deduce que 𝑧=𝑙 y 𝑧=𝑢. La segunda cuestión se corresponde con los ingredientes que recibe COSMO. En el apartado 4.1.2 se explicaba como se adaptan los ingredientes que recibe el solver ADMM para que el número de operaciones que realice sea el mínimo estrictamente necesario, pues muchas operaciones del algoritmo ADMM no cambiaban de una iteración a la siguiente. Esto se debe a que este solver esta diseñado directamente con MPC en mente. Para un solver genérico, esto no puede cumplirse, pues este no sabe de antemano si será ejecutado más de una vez para resolver problemas sucesivos en masa, como el caso del MPC, o simplemente debe resolver un único problema cuadrático. Por eso, los ingredientes que recibe COSMO son los parámetros del problema QP directamente, tal y como se manifestaban en el apartado 3.1.3. Además, como se pretende resolver un problema QP para cada tiempo de muestreo, es necesario actualizar las variables que dependen de los parámetros del sistema o del controlador, en este caso, 𝑞 y 𝑏. Para ello, pueden aprovecharse directamente las herramientas de las que dispone COSMO, llamando a la función update en cada iteración de la simulación. 5.2 Simulaciones frente a quadprog Antes de pasar a los experimentos finales, para cersiorarse del correcto funcionamiento del solver programado, se establece una simulación en la que se controla el mismo sistema a través de dos solvers diferentes. Por un lado, el solver desarrollado en el proyecto y por el otro, el solver nativo de Matlab para el problema QP, quadprog. Se realizan dos simulaciones, cada una frente a un sistema diferente. Como esta comparativa es un trámite para validar el funcionamiento del solver y no su eficiencia, no cabe entrar en mucho detalle. Simplemente, se espera que el sistema controlado converga hacia la referencia de una forma parecida a quadprog. Los sistemas expuestos se explican en profundidad en los siguientes apartados. Para introducirlos someramente, el primer sistema está compuesto por cuatro tanques de agua, cada uno representado por un estado y controlados mediante un par de grifos, dos señales de control. Es por tanto un sistema de pequeñas dimensiones y no se espera que alcance las restricciones activas, ni en estados ni en señales de control. El segundo sistema está formado por un conjunto de masas conectadas por muelles, cada una de ellas representada por un par de estados. Se actúa sobre las masas externas de la cadena. En este caso sí pueden darse dimensiones de mayor tamaño, además de que se espera que se alcancen las restricciones activas. Se realiza una simulación que sigue el esquema de Tabla 1. En primer lugar, para el sistema de los tanques de agua, los parámetros respectivos al MPC son (85), los parámetros respectivos al ADMM son (86) y la referencia (87). Los valores respectivos al modelo se detallan en el Capítulo 5, donde se explican los sistemas simulados en profundidad. En el Capítulo 5 también se explica el cálculo de 𝑃. Lás dos primeras gráficas (Ilustración 4, Ilustración 5) mostradas se corresponden a un ejemplo de estado y señal de control, respectivamente, a lo largo de la simulación. Como solo se pretende comprobar que el problema de regulación se resuelve satisfactoriamente, es decir, el control alcanza la referencia; no se muestran todos los estados y señales de control sino solo un ejemplo de cada uno, por claridad. Nótese también que como la simulación se hace a través de un modelo lineal, y este surge de linealización de las ecuaciones diferenciales en torno a un punto de operación, los ejes de ordenadas de las gráficas representan los valores incrementales de dicha variable en torno al punto de operación. La tercera de las gráficas, Ilustración 6, representa la diferencia entre la 𝑧 óptima obtenida por el solver quadprog y el solver ADMM desarrollado en el proyecto 𝑧𝑖𝑑𝑖𝑓𝑓=‖𝑧𝑖∗ 𝑄𝑈𝐴𝐷𝑃𝑅𝑂𝐺 − 𝑧𝑖∗ 𝐴𝐷𝑀𝑀 ‖∞ (81) Recuérdese que 𝑧 óptima es resultado entregado por el solver como solución del problema QP, notado como 𝑧∗. Para cada instante, se resta el vector 𝑧 obtenido de quadprog y del solver ADMM. En lugar de representar alguna 37 de las componentes individuales del vector diferencia, se representa la norma infinita, es decir, la mayor de las diferencias en el vector de la diferencia. Ilustración 4. Estado del primer tanque de agua Ilustración 5. Señal de control 1 (Grifo 1) Resultados númericos 38 Ilustración 6. Diferencia entre la 𝒛∗calculada por quadprog y el solver ADMM propio para el sistema de los tanques Para la simulación del sistema de los muelles, los parámetros son (93)-(95). Las líneas negras de las Ilustración 7 e Ilustración 8 representan la restricción superior de los respectivos estado y señal de control. La restricción inferior no se representa por claridad, pues no está ni siquiera cerca de ser alcanzada. Ilustración 7. Estado de la segunda masa 39 Ilustración 8. Señal de control 2 (Fuerza 2) Ilustración 9. Diferencia entre la 𝒛∗calculada por quadprog y el solver ADMM propio para el sistema de los muelles En efecto, puede observarse que la trayectoria seguida en todos los ejemplos mostrados por el resultado del solver ADMM desarrollado a lo largo del proyecto es prácticamente idéntica a la del solver propio de Matlab, software comercial (y, por tanto, se asume que cumple unos estándares de calidad). Resultados númericos 40 5.3 Sistema de los tanques de agua Para el primer experimento con el solver, se ha seleccionado el sistema descrito en [23]. Es un sistema compuesto por cuatro tanques de agua, de los cuales los dos superiores son llenados a través de dos bombas que pueden accionarse y suponen las señales de control del sistema. Cada uno de los tanques superiores desagua a uno de los tanques inferiores. Los estados de este sistema representan la altura que alcanza el agua en cada uno de los tanques. Se ha seleccionado este sistema para comenzar pues sus dimensiones son bastante pequeñas y es muy difícil que se alcancen las restricciones activas. Los tanques inferiores, además del desagüe, son llenados también por las válvulas de los tanques superiores, a través de un parámetro que divide en las válvulas la cantidad de caudal que reciben los tanques superiores y cuanto los inferiores. Aun así, para este problema, ese se considera un parámetro dado y viene reflejado directamente en el modelo del sistema. Para este sistema se parte directamente del modelo en espacio de estados, linealizado a partir de las ecuaciones diferenciales en torno a un punto de operación 𝑋𝑜=(0.7292 0.8102 0.6594 0.9408)𝑇,⁡ 𝑈𝑜=(1.9480 2)𝑇,𝑌𝑜=(0.7292 0.8102)𝑇 (82) que viene caracterizado por las siguientes matrices 𝐴=(0.9434 0 0.0421 0 0 0.9382 0 0.0336 0 0 0.9579 0 0 0 0 0.9664);𝐵=(0.0139 0 0 0.0185 0 0.0278 0.0324 0 ) (83) Las restricciones dadas son Ilustración 10. Diagrama esquemático del sistema de los cuatro tanques. Figura 2 en [23] 41 𝑥=(0.4708 0.3898 0.5406 0.2592)𝑇 𝑥=(−0.5292 −0.6102 −0.4594 −0.7408)𝑇 𝑢=(1.052 1)𝑇 𝑢=(−1.948 −2)𝑇 (84) Sobre este modelo se realizan las siguientes simulaciones. 5.3.1 Simulación con 𝑹=𝑰𝟐 La simualción se realiza a través de un fichero que implemente el código expuesto en Tabla 1. Por cuestiones de claridad, en lo respectivo a las gráficas que muestran estados y señales de control convergiendo a la referencia, solo se muestra un estado y una señal de control, pues estos resultados tienen poco interés más allá de demostrar el correcto funcionamiento del sistema. Los resultados numéricos interesantes de los que deben sacarse las conclusiones son los respectivos a las gráficas de tiempos e iteraciones. Los parámetros respectivos al MPC son 𝑁=10;⁡𝑄=(10 000 010 0 0 0 0 10 0 00010);⁡𝑅=(1 0 0 1);𝑃=𝑑𝑙𝑞𝑟(𝐴,𝐵,𝑄,𝑅) (85) Los parámetros respectivos al ADMM son 𝜌=15;𝜀𝑝=𝜀𝑑=10−5;⁡𝑘𝑚𝑎𝑥=500 (86) y la referencia 𝑢𝑟=(0.2 −0.1);𝑥𝑟=(𝐼−𝐴)−1·𝐵·𝑢𝑟 (87) El cálculo de 𝑃 se basa en las ecuaciones de LQR para sistemas discretos. Esta matriz representa la solución de la ecuación de Riccati para la función de coste del problema LQR. Los resultados de este experimento pueden observarse en las siguientes gráficas. La línea roja representa la referencia. Ilustración 11. Segundo estado para el sistema de los 4 tanques Resultados númericos 48 Ilustración 22. Iteraciones del solver para el sistema de las masas y los muelles sin restricciones (N=100) Ilustración 23. Tiempos de ejecución para el sistema de las masas y los muelles sin restricciones (N=100) 5.4.3 Simulación con restricciones activas Este experimento es el más interesante de todos los realizados. Esto se debe a que, como ya se ha comentado anteriormente en este documento, las restricciones suponen el elemento diferenciador entre MPC y otros sistemas de control similares como LQR. Para esta simulación, los parámetros son ⁡𝑁=10;𝑃=𝑑𝑙𝑞𝑟(𝐴,𝐵,𝑄,𝑅),⁡ 𝑄= ( 15 0 0 0 0 0 015 0 0 0 0 0 0 15 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ) ;𝑅=(0.1 0 0 0.1), (93) 49 𝜌=15;𝜀𝑝=𝜀𝑑=10−5;𝑘𝑚𝑎𝑥=1000 (94) 𝑢𝑟=(0.5 0.5)𝑇;𝑥𝑟=(0.25 0.25 0.25 1 1 1)𝑇 (95) Los resultados se muestran a continuación. En las iteraciones de Matlab y Julia puede observarse como en un par determinado de instantes de muestreo las iteraciones se desploman hasta alcanzar un valor por debajo de 0. Esto indica que para ese par de iteraciones el solver alcanzó el nivel máximo de iteraciones. A pesar de ello, el controlador puede mantener un funcionamiento adecuado y alcanzar la referencia. Por ese mismo motivo, se han mantenido los resultados del experimento en lugar de reptirlos aumentando el valor máximo admisible para las iteraciones. Las líneas negras de las siguientes gráficas representan las restricciones superiores. Las restricciones inferiores no se representan pues en ningún momento del experimento son alcanzadas. Ilustración 24. Segundo estado para el sistema de las masas y los muelles con restricciones Ilustración 25. Segunda señal de control para el sistema de las masas y los mulles con restricciones Resultados númericos 50 Ilustración 26. Iteraciones del solver para el sistema de las masas y los muelles con restricciones Ilustración 27. Tiempos de ejecución para el sistema de las masas y los muelles con restricciones 5.5 Datos estadísticos de los tiempos e iteraciones Como conclusión a los experimentos realizados, se incluye una tabla (Tabla 9) en la que se recogen la media, mediana, máximo y mínimo de todos los tiempos e iteraciones para cada uno de los solvers en cada uno de los experimentos. Los tiempos se encuentran en milisegundos a no ser que se especifíque lo contrario. Los números mostrados en rojos representan números afectados porque se alcanzó 𝑘𝑚𝑎𝑥. Una iteración marcada en rojo indica que esa es la iteración mayor de la simulación sin tener en cuenta en las que se alcanzó 𝑘𝑚𝑎𝑥. Un tiempo marcado en rojo indica que ese tiempo se obtuvo en una iteración en la que se alcanzó 𝑘𝑚𝑎𝑥. 51 Simulación (Tanques de agua) 𝑅=𝐼2 𝑅=10·𝐼2 Solver Matlab Julia COSMO Matlab Julia COSMO Iteraciones Media 80.2 80.2 15.4 19 19 13.9 Mediana 78 78 12 19 19 12 Máximo 104 104 32 19 19 25 Mínimo 78 78 5 19 19 10 Tiempos Media 25.4 34.4 6.9 6 9.1 7.7 Mediana 24.8 33.7 6.8 6 9 7.7 Máximo 34 46.2 12.1 7.5 10.8 10.5 Mínimo 24.1 30.8 3.4 5.9 8.2 5.7 Simulación (Masas y muelles ) Sin restricciones activas y 𝑁=10 Sin restricciones activas y 𝑁=100 Con restricciones activas Solver Matlab Julia COSMO Matlab Julia COSMO Matlab Julia COSMO Iteraciones Media 431.7 431.7 8.3 431.2 431.2 11.36 534.2 534.2 21.8 Mediana 425 425 7.5 425 425 10 526 526 10 Máximo 633 633 52 632 632 55 656 656 151 Mínimo 414 414 5 425 425 5 517 517 5 Tiempos Media 337.6 484.8 9.3 3.38s 5.11s 11.1 448.8 627.4 10.2 Mediana 335.5 482.2 9.5 3.37s 5.08 11.4 428.1 598.5 7.6 Máximo 417.5 596.7 11.3 4.19s 6.29s 15.3 814.3 1139 40.3 Mínimo 326.1 468.8 7.4 3.34s 5.06s 7.5 419.5 584.2 5.8 Tabla 9. Datos estadísticos de los tiempos e iteraciones 53 6 CONCLUSIONES Los experimentos suponen la conclusión del proyecto. Con los resultados numéricos expuestos, la etapa final de esta memoria es el comentario y la reflexión que suscitan dichos resultados. Aun así, previamente, también cabe una reflexión sobre todo el recorrido que ha supuesto desarrollar este proyecto y como ha influido en mi formación como ingeniero. El punto principal a tratar en este proyecto, sobre todo desde el punto de vista formativo, es MPC. Uno de los ejes vertebradores del currículo del Grado en Ingeniería Electrónica, Robótica y Mecatrónica es la automatización y los sistemas de control. Aun así, se trata de una disciplina enorme con muchísimas vertientes y aplicaciones, por lo que el grueso de los estudios está compuesto por las cuestiones más fundamentales de la materia. Queda entonces poco espacio para incluir conocimientos sobre técnicas más avanzadas y punteras en el sector. Por eso mismo, valoro mucho el haberme podido acercar a una técnica de control bastante más avanzada e interesante como es el control predictivo, y en concreto, MPC. Aunque también se trata de un conocimiento muy amplio, imposible de abarcar al completo en un proyecto de este calibre, he podido cimentar muy bien las bases del mismo y adquirir las herramientas necesarias para continuar formándome en un campo tan prometedor. Además, me gustaría resaltar que, al tener que adaptar MPC a ADMM, no solo he aprendido las superficialidades de la técnica, sino que me ha permitido trabajar más en profundidad sobre algunos puntos. El punto de intersección entre MPC y ADMM, la optimización matemática, es el otro principal campo sobre el que he podido formarme. Se trata de otro campo en el que también la formación estándar abarca solo los puntos más básicos y solo a través de trabajo individual se adquieren conocimientos más avanzados. He podido familiarizarme con los conceptos que vertebran la disciplina y he podido poner en práctica los conocimientos matemáticos que he adquirido a lo largo del grado. Dentro de la optimización, además de conocer de manera generalista las ideas clave, he podido enfrentarme a un problema y solución reales en la forma del algoritmo ADMM. Estos puntos constituyen de manera general los conocimientos teóricos que he tenido la oportunidad de adquirir. Por otro lado, se encuentra el trabajo práctico desarrollado. La computación y los lenguajes de programación han sido uno de los mayores estímulos que he tenido durante estos años. Este proyecto supone el culmen de la puesta a punto de mis conocimientos sobre este campo. Además de enfrentarme a un problema de mayor envergadura en Matlab, he tenido la ocasión de aprender las bases de un lenguaje muy puntero y prometedor, Julia. Además, he desarrollado una herramienta completa, compuesta de un conjunto de códigos estructurados y diseñados para obtener la mayor eficiencia y modularidad posible. También he realizado diversas pruebas de la misma, con resultados muy satisfactorios. El solver programado, más allá de las aptitudes que tiene comparadas a las de otros solvers, es completamente funcional y ha quedado validado a través de los experimentos. Aun así, es importante señalar que los resultados de eficiencia no se corresponden con lo esperado. A pesar de que los resultados obtenidos en simulación de sistemas suficientemente sencillos sean buenos, los resultados frente a sistemas mayores, dentro de que son aceptables, no pueden compararse con los de un solver extremadamente competitivo como es COSMO. Además, para sistemas aun mayores, los tiempos empeoran bastante, comparados con el buen escalado que mantiene COSMO, con lo que el escalado del solver no es un punto a destacar del mismo. Aun así, el principal fracaso ha sido la eficacia esperada de Julia. Se pretendía que la implementación en Julia supusiese una mejora sustancial en los tiempos de ejecución del solver, en base a las características particulares del sistema. Sin embargo, no solo ha sido así, sino que se ni siquiera se han mejorado los tiempos en Matlab. Sospecho que la principal causa de este inconveniente esté relacionada con el hecho de realizar la implementación en Julia a través de una transcripción de Matlab. Conocer en mayor profundidad los entresijos del lenguaje y dedicar más tiempo a estudiar la funcionalidad de las características que lo hacen tan veloz, podría haber dado lugar a una implementación que explotase estás ventajas en mayor profundidad. En Conclusiones 54 concreto, la funcionalidad de multiple dispatch es un elemento diferenciador de Julia que colabora mucho a hacerlo tan rápido y puede ser un de los puntos a explotar para mejorar el rendimiento. Aun así, considero un acierto el procedimiento seguido para la realización de este proyecto, quizá no tanto desde el aspecto tecnológico, sino desde el formativo, pues esta línea de desarrollo me ha permitido profundizar más en los conocimientos sobre MPC, ADMM, álgebra sparsa y Matlab. A pesar de todo, los puntos débiles de los resultados obtenidos pueden presentarse como una oportunidad de futuro. Mejorar los códigos que constituyen este proyecto suponen tanto un reto teórico, en cuanto conocer y aprender mejor sobre el interior de un lenguaje de programación y como estos funcionan, de cara a sacarle el mayor partido posible al solver; como un reto práctico, a la hora de poner estas mejoras en funcionamiento y validarlas, obtenido mejores resultados en los tiempos de ejecución del solver. El otro camino que se abre paso de manera evidente al finalizar este proyecto es algo que ya se ha mencionado en varias ocasiones. El enfoque con el que se perfila este proyecto está estrechamente relacionado con el trabajo de Pablo Krupa, uno de los profesores encargados de la supervisión de este proyecto. Krupa ha estado contribuyendo activamente a la investigación en torno a MPC y una de las líneas de investigación que se encuentra siguiendo es la de adaptar esta técnica tan computacionalmente costosa a los dispositivos más modestos y con menores recursos posibles. Esto es una oportunidad importante de fomentar el uso de MPC en la industria. De aquí surge la preocupación que se da en este proyecto a la eficiencia. Es por eso que el acercamiento a la utilidad real del desarrollo del solver es su implementación en dispositivos humildes, como en un PLC o en sistemas embebidos. Como trabajo futuro queda la adaptación del solver aquí desarrollado a un sistema real controlado a través de un dispositivo de escasos recursos y la validación del mismo en un entorno real. 55 ANEXO A: CÓDIGOS DE LA IMPLEMENTACIÓN DEL SOLVER EN MATLAB Generatriz (gen_ADMM_MPC.m) function out=gen_ADMM_MPC(model,rho) 1 2 3 %%% Dimensiones 4 n=model.n; 5 m=model.m; 6 N=model.N; 7 dimz=N*(n+m); 8 9 %%% Cálculo de las restricciones 10 indx=0; 11 LB=zeros(dimz,1); 12 UB=zeros(dimz,1); 13 while indx<dimz 14 15 for i=1:m 16 LB(indx+i)=model.LBu(i); 17 UB(indx+i)=model.UBu(i); 18 end 19 indx=indx+m; 20 21 for i=1:n 22 LB(indx+i)=model.LBx(i); 23 UB(indx+i)=model.UBx(i); 24 end 25 indx=indx+n; 26 27 end 28 29 %%% Cálculo de Aeq 30 i=0; 31 j=0; 32 t=0; 33 I=-eye(n); 34 A=zeros(N*n,dimz); 35 while i<(N*n) 36 37 for di=1:n 38 for jA=1:t 39 A(i+di,j+jA)=model.A(di,jA); 40 end 41 for jB=(t+1):(t+m) 42 A(i+di,j+jB)=model.B(di,jB-t); 43 end 44 for jI=(t+m+1):(t+m+n) 45 56 A(i+di,j+jI)=I(di,jI-(t+m)); 46 end 47 end 48 49 i=i+n; 50 j=j+t+m; 51 t=n; 52 53 end 54 55 %%% Cálculo de Hz y W 56 H=[]; 57 for i=1:N-1 58 H=blkdiag(H,model.R,model.Q); 59 end 60 H=blkdiag(H,model.R,model.P); 61 62 Hz=H+rho*eye(size(H,1)); 63 Hi=inv(Hz); 64 AHi=A*Hi; 65 HiA=Hi*A'; 66 W=AHi*A'; 67 68 %%% Transformación a sparsa 69 Hi_csr=full2CSR(Hi); 70 AHi_csr=full2CSR(AHi); 71 HiA_csr=full2CSR(HiA); 72 [W_L,W_Di]=sparsa_decomp(W); 73 74 %%% Salidas 75 out.LB=LB; 76 out.UB=UB; 77 78 out.Hi_csr=Hi_csr; 79 out.AHi_csr=AHi_csr; 80 out.HiA_csr=HiA_csr; 81 out.W_L=W_L; 82 out.W_Di=W_Di; 83 84 85 end 86 57 Solver (solver_ADMM_MPC.m) function out=solver_ADMM_MPC(model,x_act,p,var) 1 2 3 %%% Cálculo de q y b 4 n=model.n; 5 m=model.m; 6 N=model.N; 7 dimz=N*(n+m); 8 9 qR = -model.R*model.u_r; 10 qQ = -model.Q*model.x_r; 11 qP = -model.P*model.x_r; 12 q_aux=[qR;qQ]; 13 q=zeros(dimz,1); 14 for i=1:(n+m):(dimz-n-m) 15 q(i:i+n+m-1)=q_aux; 16 end 17 q_aux=[qR;qP]; 18 q(dimz-n-m+1:dimz)=q_aux; 19 20 b_aux=-model.A*x_act; 21 b=zeros(N*n,1); 22 b(1:n)=b_aux; 23 24 25 %%% Inicialización del sistema 26 v_k=zeros(dimz,1); 27 v_k1=zeros(dimz,1); 28 lmbda_k=zeros(dimz,1); 29 finish=false; 30 k=1; 31 32 ri=1/p.rho; 33 while ~finish 34 35 %%% Cálculo de z (sparsa) 36 qz_k=q+lmbda_k-p.rho*v_k; 37 38 b_hat = -b - sparsa_prod(var.AHi_csr,qz_k); 39 mu_k=sparsa_ecsys(var.W_L,var.W_Di,b_hat); 40 41 num1=sparsa_prod(var.HiA_csr,mu_k); 42 num2=sparsa_prod(var.Hi_csr,qz_k); 43 z_k1 = -num1 - num2; 44 45 %%% Cálculo de v 46 qv_k=lmbda_k+p.rho*z_k1; 47 for i=1:dimz 48 v_k1(i)=max(min(qv_k(i)*ri,var.UB(i)),var.LB(i)); 49 end 50 51 %%% Cálculo de lambda 52 lmbda_k1=lmbda_k+p.rho*(z_k1-v_k1); 53 54 %%% Conidción de salida 55 64 57 m_MPC.N=10; 58 m_MPC.Q=Q; 59 m_MPC.R=R; 60 [~,m_MPC.P]=dlqr(A,B,Q,R); 61 62 %%% Parámetros ADMM 63 p_ADMM.rho=15; 64 p_ADMM.k_max=k_max; 65 p_ADMM.eps_p=1e-5; 66 p_ADMM.eps_d=1e-5; 67 68 %%% Simulación 69 eTic=tic; 70 for i=1:n_evals 71 sTic=tic; 72 hist=sim_system(m_MPC,p_ADMM,samples); 73 t_sim(i)=toc(sTic); 74 t_j=t_j+1; 75 end 76 t_ADMM=toc(eTic); 77 t_sim_ADMM=sum(t_sim); 78 t_sol_ADMM=sum(sum(t_sol)); 79 80 t_sol=1000*t_sol; 81 t_mean_ADMM=zeros(samples,1); 82 for i=1:samples 83 t_mean_ADMM(i)=sum(t_sol(i,:))/double(n_evals); 84 end 85 hist.time=t_mean_ADMM; 86 87 %%% Gráficas 88 if plots 89 if springs 90 run('plot_6Spring.m'); 91 else 92 run('plot_4Tank.m'); 93 end 94 end 95 96 save('results.mat','hist'); 97 65 ANEXO B: CÓDIGOS DE LA IMPLEMENTACIÓN DEL SOLVER EN JULIA Generación de espacio de trabajo en Julia (Generate_Julia_Workspace.m) 1 2 springs=0; 3 4 %%% Datos 5 if springs 6 load('6Spring.mat') 7 s=sys; 8 u_r=[0.5;0.5]; 9 samples=50; 10 else 11 load('4Tank.mat') 12 s=sys; 13 u_r=[0.2;-0.1]; 14 samples=150; 15 end 16 17 MPC.A=s.A; 18 MPC.B=s.B; 19 MPC.LBx=s.LBx; 20 MPC.UBx=s.UBx; 21 MPC.LBu=s.LBu; 22 MPC.UBu=s.UBu; 23 MPC.u_r=u_r; 24 25 %%% Cálculo del estado de referencia 26 if springs 27 x_r=0.25*[1;1;1;0;0;0]; 28 else 29 x_r=(eye(size(s.A))-s.A)\s.B*u_r; 30 end 31 MPC.x_r=x_r; 32 33 %%% Parámetros MPC 34 if springs 35 Q=diag([15,15,15,1,1,1]); 36 R=0.1*eye(2); 37 else 38 Q=10*eye(4); 39 R=10*eye(2); 40 end 41 42 MPC.N=10; 43 MPC.Q=Q; 44 MPC.R=R; 45 [~,MPC.P]=dlqr(MPC.A,MPC.B,Q,R); 46 47 %%% Parámetros ADMM 48 66 ADMM.rho=15; 49 ADMM.k_max=500; 50 ADMM.eps_p=1e-5; 51 ADMM.eps_d=1e-5; 52 53 %%% Parámetros simulación 54 [n,m]=size(MPC.B); 55 MPC.n=n; 56 MPC.m=m; 57 mat_COSMO=gen_COSMO_MPC(MPC,zeros(n,1)); 58 mat_ADMM=gen_ADMM_MPC(MPC,ADMM.rho); 59 60 %%% Generación de fichero 61 data.model=MPC; 62 data.param=ADMM; 63 data.samples=samples; 64 data.matrixes=mat_COSMO; 65 save('WSJulia_COSMO.mat','data'); 66 67 data.matrixes=mat_ADMM; 68 save('WSJulia_ADMM.mat','data'); 69 NOTA: La función generatriz, según lo explicado en el apartado 4.2.1, es utilizada directamente desde Matlab debido a la modularidad de la solución propuesta. Por tanto, no se expone en este anexo, sino que se remite al lector a la función correspondiente del Anexo A, pues es idéntica a la de la implementación de Julia. Lo mismo ocurre para la función que realiza la descomposición LDL sparsa y para los conversores CSR y CSC. 67 Generatriz para COSMO (gen_COSMO_MPC.m) function out=gen_COSMO_MPC(model,x_act) 1 2 3 %%% Dimensiones 4 n=model.n; 5 m=model.m; 6 N=model.N; 7 dimz=N*(n+m); 8 9 %%% Cálculo de las restricciones 10 indx=0; 11 LB=zeros(dimz,1); 12 UB=zeros(dimz,1); 13 while indx<dimz 14 15 for i=1:m 16 LB(indx+i)=model.LBu(i); 17 UB(indx+i)=model.UBu(i); 18 end 19 indx=indx+m; 20 21 for i=1:n 22 LB(indx+i)=model.LBx(i); 23 UB(indx+i)=model.UBx(i); 24 end 25 indx=indx+n; 26 27 end 28 29 %%% Cálculo de Aeq 30 i=0; 31 j=0; 32 t=0; 33 I=-eye(n); 34 A=zeros(N*n,dimz); 35 while i<(N*n) 36 37 for di=1:n 38 for jA=1:t 39 A(i+di,j+jA)=model.A(di,jA); 40 end 41 for jB=(t+1):(t+m) 42 A(i+di,j+jB)=model.B(di,jB-t); 43 end 44 for jI=(t+m+1):(t+m+n) 45 A(i+di,j+jI)=I(di,jI-(t+m)); 46 end 47 end 48 49 i=i+n; 50 j=j+t+m; 51 t=n; 52 53 end 54 55 68 %%% Cálculo de H 56 H=[]; 57 for i=1:N-1 58 H=blkdiag(H,model.R,model.Q); 59 end 60 H=blkdiag(H,model.R,model.P); 61 62 %%% Cálculo de q y b 63 qR = -model.R*model.u_r; 64 qQ = -model.Q*model.x_r; 65 qP = -model.P*model.x_r; 66 q_aux=[qR;qQ]; 67 q=zeros(dimz,1); 68 for i=1:(n+m):(dimz-n-m) 69 q(i:i+n+m-1)=q_aux; 70 end 71 q_aux=[qR;qP]; 72 q(dimz-n-m+1:dimz)=q_aux; 73 74 b_aux=-model.A*x_act; 75 b=zeros(N*n,1); 76 b(1:n)=b_aux; 77 78 %%% Salidas 79 out.LB=LB; 80 out.UB=UB; 81 82 out.H=H; 83 out.A=A; 84 out.q=q; 85 out.b=b; 86 87 end 88 69 Solver ADMM (solver_ADMM.jl) using LinearAlgebra 1 include("solver_ADMM_aux.jl") 2 3 4 function solver_ADMM_MPC(model,x_act,p,var) 5 6 ### update_qb 7 n=convert(Int64,model["n"]); 8 m=convert(Int64,model["m"]); 9 N=convert(Int64,model["N"]); 10 dimz=N*(n+m); 11 12 qR = -model["R"]*model["u_r"]; 13 qQ = -model["Q"]*model["x_r"]; 14 qP = -model["P"]*model["x_r"]; 15 q_aux=[qR;qQ]; 16 q=zeros(dimz); 17 for i in 1:(n+m):(dimz-n-m) 18 q[i:i+n+m-1]=q_aux; 19 end 20 q_aux=[qR;qP]; 21 q[dimz-n-m+1:dimz]=q_aux; 22 23 b_aux = -model["A"]*x_act; 24 b=zeros(N*n); 25 b[1:n]=b_aux; 26 27 ### Inicialización del sistema 28 v_k=zeros(dimz); 29 v_k1=zeros(dimz); 30 lmbda_k=zeros(dimz); 31 finish=false; 32 k=1; 33 34 ri=1/p["rho"]; 35 while !finish 36 37 ### Cálculo de z (sparsa) 38 qz_k=q+lmbda_k-p["rho"]*v_k; 39 40 b_hat = -b - sparsa_prod(var["AHi_csr"],qz_k); 41 mu_k=sparsa_ecsys(var["W_L"],var["W_Di"],b_hat); 42 43 num1=sparsa_prod(var["HiA_csr"],mu_k); 44 num2=sparsa_prod(var["Hi_csr"],qz_k); 45 global z_k1 = -num1 - num2; 46 47 ### Cálculo de v 48 qv_k=lmbda_k+p["rho"]*z_k1; 49 for i=1:dimz 50 UBmin=var["UB"][i,1]; 51 LBmax=var["LB"][i,1]; 52 v_k1[i]=max(min(qv_k[i]*ri,UBmin),LBmax); 53 end 54 70 55 ### Cálculo de lambda 56 lmbda_k1=lmbda_k+p["rho"]*(z_k1-v_k1); 57 58 ### Condición de salida 59 global k_bool=p["k_max"]<k; 60 fc_p=norm(z_k1-v_k1,Inf)<p["eps_p"]; 61 fc_d=norm(v_k1-v_k,Inf)<p["eps_d"]; 62 finish = (fc_p && fc_d) || k_bool; 63 64 ### Actualización de variables 65 k=k+1; 66 v_k=copy(v_k1); 67 lmbda_k=copy(lmbda_k1); 68 69 end 70 71 ### Salidas 72 if k_bool 73 k=0; 74 end 75 76 sol=z_k1; 77 out=Dict("u"=> sol[1:m], 78 "z"=> sol, 79 "k"=>k-1); 80 return out 81 82 end 83 71 Funciones auxiliares para el solver ADMM (solver_ADMM_aux.jl) using LinearAlgebra 1 2 3 function sparsa_prod(M,v) 4 5 x=zeros(convert(Int64,M["nrow"])); 6 for i=1:convert(Int64,M["nrow"]) 7 r_start=convert(Int64,M["row"][i]); 8 r_end=convert(Int64,M["row"][i+1]); 9 10 for j=r_start:r_end-1 11 x[i]=x[i]+M["val"][j]*v[convert(Int64,M["col"][j])]; 12 end 13 end 14 return x 15 16 end 17 18 19 function sparsa_ecsys(L,Di,b) 20 21 # Algoritmo QDLDL 22 x=b; 23 n=size(Di,1); 24 25 for i=1:n # Forward substitution 26 c_start=convert(Int64,L["col"][i]); 27 c_end=convert(Int64,L["col"][i+1]); 28 29 for j=c_start:c_end-1 30 Lrow=convert(Int64,L["row"][j]) 31 x[Lrow]=x[Lrow]-L["val"][j]*x[i]; 32 end 33 end 34 35 for i=1:n 36 x[i]=x[i]*Di[i]; 37 end 38 39 for i=n:-1:1 # Backwards substitution 40 c_start=convert(Int64,L["col"][i]); 41 c_end=convert(Int64,L["col"][i+1]); 42 43 for j=c_start:c_end-1 44 Lrow=convert(Int64,L["row"][j]) 45 x[i]=x[i]-L["val"][j]*x[Lrow]; 46 end 47 end 48 49 return x 50 51 end 52 72 Simulación del sistema con el solver ADMM (solver_sim_ADMM.jl) using MAT, LinearAlgebra 1 include("solver_ADMM.jl") 2 3 4 ### Lectura de los datos de Matlab 5 data=matread("WSJulia_ADMM.mat") 6 7 model=data["data"]["model"]; 8 n=convert(Int64,model["n"]); 9 m=convert(Int64,model["m"]); 10 N=convert(Int64,model["N"]); 11 12 param=data["data"]["param"]; 13 14 matrixes=data["data"]["matrixes"]; 15 samples=convert(Int64,data["data"]["samples"]); 16 17 18 ### Inicialización del sistema 19 dimz=N*(n+m); 20 hist=Dict("x"=>zeros(n,samples), 21 "u"=>zeros(m,samples), 22 "z"=>zeros(dimz,samples), 23 "k"=>zeros(1,samples), 24 "x_r"=>model["x_r"], 25 "u_r"=>model["u_r"], 26 "t_sol"=>zeros(samples)); 27 x_k=zeros(n,1); 28 29 30 t_sol=zeros(samples); 31 for i in 1:samples 32 33 ### Resolución del problema cuadrático 34 t_aux=time(); 35 sol=solver_ADMM_MPC(model,x_k,param,matrixes); 36 t_sol[i]=time()-t_aux; 37 38 39 ### Simulación del sistema 40 x_k1=model["A"]*x_k+model["B"]*sol["u"]; 41 42 43 ### Actualización de variables y guardado de datos 44 hist["x"][:,i]=x_k; 45 hist["u"][:,i]=sol["u"]; 46 hist["z"][:,i]=sol["z"]; 47 hist["k"][i]=sol["k"]; 48 global x_k=x_k1; 49 50 end 51 52 53 ### Escritura de resultados para Matlab 54 73 hist["t_sol"]=t_sol; 55 matwrite("WSJulia_ADMM_results.mat",hist); 56 Funciones auxiliares para el solver COSMO (solver_COSMO_aux.jl) using LinearAlgebra 1 2 3 function update_qb(model,x_act) 4 n=convert(Int64,model["n"]); 5 m=convert(Int64,model["m"]); 6 N=convert(Int64,model["N"]); 7 dimz=N*(n+m); 8 9 ### Cálculo de q y b (provisional) 10 qR = -model["R"]*model["u_r"]; 11 qQ = -model["Q"]*model["x_r"]; 12 qP = -model["P"]*model["x_r"]; 13 q_aux=[qR;qQ]; 14 q=zeros(dimz); 15 for i in 1:(n+m):(dimz-n-m) 16 q[i:i+n+m-1]=q_aux; 17 end 18 q_aux=[qR;qP]; 19 q[dimz-n-m+1:dimz]=q_aux; 20 21 b_aux = -model["A"]*x_act; 22 b=zeros(N*n); 23 b[1:n]=b_aux; 24 25 return q, b 26 end 27 80 [15] M. Garstka, M. Cannon y P. Goulart, «COSMO: A conic operator splitting method for large convex problems,» de European Control Conference, pp. 1951-1956, Nápoles, Italia, 2019. [16] S. Boyd et al., «Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,» Foundations and Trends© in Machine Learning, vol. 3, nº 1, pp. 1-122, 2011. [17] H. Kwakernaak y R. Sivan, Linear Optimal Control Systems, John Wiley & Sons Inc., 1972. [18] S. Boyd y L. Vandenberghe, «Convex Optimization,» Cambridge University Press, 2004. [19] A. Buluç, J. Fineman, M. Frigo, J. Gilbert y C. Leiserson, «Parallel Sparse Matrix-Vector and MatrixTranspose-Vector Multiplication Using Compressed Sparse Blocks,» de ACM Symposium on Parallelism in Algorithms and Architectures, pp. 233-244, Calgary, Alberta, Canada, 2009. [20] B. Stellato, G. Banjac, P. Goulart, A. Bemporad y S. Boyd, «OSQP: an operator splitting solver for quadratic programs,» Mathematical Programming Computation, vol. 12, nº 4, pp. 637-672, 2020. [21] JuliaLang, «Julia 1.6 Documentation—Scope of Variables,» [En línea]. Available: https://docs.julialang.org/en/v1/manual/variables-and-scoping/. [Último acceso: Junio 2021]. [22] University of Oxford Control Group, «COSMO.jl - A quadratic objective conic solver implemented in pure Julia,» [En línea]. Available: https://github.com/oxfordcontrol/COSMO.jl. [Último acceso: Junio 2021]. [23] K. Johansson, «The quadruple-tank process: a multivariable laboratory process with an adjustable zero,» IEEE Transactions on Control Systems Technology, vol. 8, nº 3, pp. 456-465, 2000. [24] P. Krupa, R. Jaouani, D. Limon y T. Alamo, «A sparse ADMM-based solver for linear MPC subject to terminal quadratic constraint,» Mayo 2021. [En línea]. Available: https://arxiv.org/abs/2105.08419. [Último acceso: Junio 2021].