Procesamiento masivamente paralelo en control predictivo basado en datos

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Procesamiento masivamente paralelo en control predictivo basado en datos Autor: Alfonso Daniel Carnerero Panduro Tutor: Daniel Rodríguez Ramírez Dep. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Procesamiento masivamente paralelo en control predictivo basado en datos Autor: Alfonso Daniel Carnerero Panduro Tutor: Daniel Rodríguez Ramírez Profesor Titular Dep. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 Trabajo Fin de Grado: Procesamiento masivamente paralelo en control predictivo basado en datos Autor: Alfonso Daniel Carnerero Panduro Tutor: Daniel Rodríguez Ramírez El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos A mis familiares, amigos y profesores. I Resumen Durante mucho tiempo, la computación paralela ha estado relegada a supercomputadores de grandes centros de investigación. Sin embargo, con la aparición de la arquitectura CUDA en las tarjetas gráficas de NVIDIA (las cuales no suponían una inversión desorbitada), empezó a estar disponible para un público mucho más amplio. Por otro lado, comprobaremos que la programación con las APIs proporcionadas por el CUDA Toolkit se plantea realmente muy similar a la que estamos acostumbrados a usar habitualmente, permitiéndonos usar diversos lenguajes populares como C o Fortran. La razón por la que planteamos el estudio de estas técnicas encuentra su sentido en el hecho de que existen gran cantidad de problemas que pueden beneficiarse del uso de este tipo de soluciones, más ahora que estamos moviendo continuamente cantidades enormes de información. Dicho esto, el objetivo es explorar las distintas posibilidades actuales para programar en la GPU usando software comercial como Matlab y, una vez hecho esto, implementar algoritmos que usen esta tecnología. En concreto, se resuelven principalmente 2 problemas: •Uso de una base de datos para predecir el estado de un sistema dinámico. • Resolución de un problema de control predictivo multivariable usando una base de datos en lugar de un modelo matemático del sistema. También comentaremos las ventajas e inconvenientes de este tipo de soluciones, haciendo especial hincapié en las mejoras de rendimiento que conseguiremos con respecto a las versiones de los algoritmos anteriores sin paralelizar. III XÍndice 6.2 Controladores predictivos basados en datos 47 6.2.1 Algoritmo de control 49 6.2.2 Tareas en la GPU 52 6.2.3 Resultados de simulación 52 Índices de desempeño 55 Factibilidad 55 Coste computacional y mejora 56 Hebras y bloques 57 6.2.4 Pruebas en la planta real 57 7 Conclusiones y posibles avances 61 Índice de Figuras 63 Índice de Tablas 65 Bibliografía 67 1 Introducción Este trabajo intenta ser una introducción al mundo de la programación paralela de manera que a partir de los ejemplos que aquí se plantean quede claro que la implementación de este tipo de técnicas no es una tarea tan complicada como ha venido siendo durante mucho tiempo y que existen cantidad de potenciales aplicaciones. El esquema aproximado que seguiremos a lo largo de este escrito será el siguiente: 1. En primer lugar, hablaremos de la historia de la arquitectura CUDA y de sus principales características. Esto incluye el flujo de información en la tarjeta gráfica, los tipos de memoria que posee, la estructura de sus hilos, etc. 2. Posteriormente, propondremos un ejemplo de programación básico para poder comparar las diferencias entre un programa en la CPU y otro en la GPU. Además, veremos la problemática asociada a este tipo de soluciones, lo cual nos impulsa a centrarnos en unos problemas determinados. 3. Más adelante, estudiaremos algunas posibilidades incluidas en el Matlab Parallel Computing Toolbox. 4. Una vez que conozcamos cómo funcionan las herramientas que vamos a utilizar, plantearemos el problema de modelado de un sistema dinámico usando únicamente una base de datos, sin ningún tipo de funciones analíticas para definir su comportamiento. Comprobaremos a su vez la bondad de este último. 5. Una vez resuelto el problema de modelado, solucionaremos el problema de control de la planta de 4 tanques con un algoritmo basado igualmente en bases de datos, haciendo un estudio detallado del comportamiento y de su rendimiento. 6. Finalmente, sacaremos unas conclusiones finales en función de los resultados obtenidos. 1 2 Unidades de procesamiento gráfico y CUDA En la página web de NVIDIA se define CUDA como: "...una arquitectura de cálculo paralelo de NVIDIA que aprovecha la gran potencia de la GPU (unidad de procesamiento gráfico) para proporcionar un incremento extraordinario del rendimiento del sistema." [1] Pero, ¿qué es una GPU? Una GPU (Graphics Processor Unit) no es más que un coprocesador encargado de aligerar la carga de la CPU (Central Processing Unit) en algunas tareas determinadas. En el caso de nuestro día a día, se encarga de presentar los gráficos por pantalla. Antes de pasar a estudiar su estructura y diferencias con la CPU podría ser interesante esbozar su evolución histórica. 2.1 Historia de la CPU y GPU Durante mucho tiempo, la manera de mejorar el rendimiento de un ordenador personal estaba completamente enfocado a la mejora de su CPU. Así, durante 30 años los fabricantes fueron siendo capaces de aumentar la velocidad de sus procesadores a través de mejoras en los relojes internos de estos. De esta manera, pasaron de funcionar a unas frecuencias en torno a 1MHz a 1 GHz. Esto suponía mejorar la velocidad 1000 veces. Pero, debido a las limitaciones que poseen los procesos de fabricación, el tamaño mínimo que transistores deben tener en un chip y problemas relacionados con la temperatura, los fabricantes e investigadores tuvieron que empezar a buscar otras posibilidades. Paralelamente, se encontraba el terreno de los supercomputadores. Estos se estaban beneficiando de las mejoras de los ordenadores personales pero, empezaron a desarrollar otra tecnología que terminó resultando de gran interés para los investigadores y fabricantes anteriormente mencionados. No se trataba de otra cosa que la de poner a trabajar varios procesadores en paralelo. De esta manera, las CPUs empezaron a a incluir varios procesadores independientes, cosa a la que ya estamos bastante acostumbrados actualmente. ¿Y dónde entra en juego la GPU? Pues bien, estas empezaron a desarrollarse sobre los años 80 como consecuencia del éxito de sistemas operativos con interfaz gráfica. En aquel entonces, algunos usuarios empezaron a demandar dispositivos que mejoraran el rendimiento del sistema a la hora de presentar los gráficos en pantalla. De igual manera que dos mercados muy distintos estuvieron trabajando con las CPUs de manera independiente, esto se repitió en el terreno de las GPUs. Por ello, en el terreno de la computación más profesional, la compañía Silicon Graphics empezó a popularizar el uso de gráficos en 3D con gran cantidad de aplicaciones. En el año 1992, liberaron su interfaz de programación y sus librerías (OpenGL) pretendiendo que se convirtiera en el estándar de creación de gráficos tridimensionales. Como no podía ser de otra manera, los consumidores empezaron a demandar esta tecnología debido sobre todo al lanzamiento de videjuegos inmersivos en primera persona como Doom oQuake. A partir de aquí, varias empresas como NVIDIA, ATI Technologies y 3dfx empezaron a desarrollar procesadores gráficos 3 4Capítulo 2. Unidades de procesamiento gráfico y CUDA más asequibles para satisfacer la demanda de estos usuarios. En paralelo al uso lúdico de las tarjetas gráficas, algunos investigadores empezaron a estudiar la posibilidad de usarlas para algo distinto a la renderización de unos determinados gráficos. Sin embargo, la única manera de lidiar con las GPUs por aquel entonces era a través de su API (Application Programming Interface). El funcionamiento era el siguiente: las aplicaciones tenían que pasar como datos a las tarjetas gráficas colores, texturas, luminosidad, etc. Lo interesante venía cuando los programadores empezaron a enviar información distinta, aunque disfrazada como la mencionada anteriormente. De esta manera, se podían realizar algunas operaciones. Sin embargo, esta programación era tan restrictiva y dependiente del modelo concreto de la GPU que no logró mucho éxito dentro del sector. Esto cambió cuando NVIDIA presentó la arquitectura CUDA (Compute Unified Device Architecture) en su modelo GeForce 8800 GTX, la cual alivia gran parte de las limitaciones de sus predecesoras en el terreno de la computación de propósito general. 2.2 Arquitectura CUDA Gracias a la cantidad de núcleos que poseen las tarjetas gráficas, existen un buen número de problemas que pueden ver mejorado su rendimiento haciendo uso de la computación paralela. A pesar de la mejora evidente en velocidad que se nos sugiere al aumentar en tal medida la cantidad de núcleos, existen alguna deficiencias como por ejemplo: • Siempre hay un tiempo implícito en cada ejecución relacionado con el paso de memoria entre la CPU y GPU. •Los núcleos CUDA deben ejecutar todos el mismo código. Es decir, no son independientes. • La sincronización entre ellos es muy restrictiva, de igual modo ocurre con los posibles intercambios de información. El flujo de procesos de una ejecución en CUDA viene representado en la figura 2.1. Figura 2.1 Representación del flujo de información en CUDA 1. 1Imagen procedente de: https://en.wikipedia.org/wiki/CUDA 2.2 Arquitectura CUDA 5 En ella vemos como es necesario mover la información varias veces durante la ejecución de un programa. Los tiempos que son necesarios para este traspaso (y que son estrictamente necesarios) son denominados overheads y son completamente decisivos a la hora de elegir usar una GPU para una determinada aplicación. Por lo general, la mejora en velocidad surge para problemas a gran escala. Por el contrario, el tiempo de computación para problemas pequeños se mantiene igual al de un programa convencional en la CPU, llegando a empeorar en algunos casos. Además, es necesario que el problema que queramos tratar sea paralelizable. Para explicar este concepto, haremos uso de varios ejemplos. • Pensemos en primer lugar en las aplicaciones para las que fueran pensadas las GPUs, es decir, procesar imágenes. Una vez tenemos la información necesario del entorno que queremos representar (luces, colores, etc) cada pixel es completamente independiente del resto por lo que es posible calcular los valores RGB que son necesarios para mostrarlo por pantalla sin necesidad de conocer el resultado del resto. No importa de ninguna manera el orden de ejecución ni existe ninguna relación de dependencia. • Planteemos ahora, en contraposición, la resolución de un sistema de ecuaciones diferenciales en tiempo discreto. En este caso no es posible resolver el sistema para un instante temporal en el que todavía no conozco el resultado del instante anterior. Es decir, debe iterarse en orden y éste no puede verse alterado. En cambio, ¿qué ocurriría si lo que quiero es resolver muchas veces ese sistema de ecuaciones en diferencias para distintos valores iniciales? En este caso vemos como cada núcleo de la GPU se encargará de realizar las operaciones convenientes de su trayectoria, completamente independiente de las del resto de núcleos. Esto significa que la GPU nos plantea dos posibilidades de aplicación: •Resolver problemas que son implícitamente paralelizables. •Resolver múltiples problemas no paralelizables. Algunos ejemplos de dispositivos que estén funcionando actualmente usando CUDA podríamos mencionar: • Un sistema llamado TechniScan Medical dedicado a la recreación de imágenes tridimensionales a partir de ultrasonidos para la prevención del cáncer de mama [2]. • Prototipos de UAVs que funcionan con NVIDIA Jetson y tienen su funcionamiento basado en algoritmos de deep-learning [3]. Como podemos ver, el campo de aplicación de CUDA es realmente diverso. La siguiente pregunta podría ser cómo está estructurada la memoria y la ejecución de programas en CUDA. De momento, nos basta con imaginarnos que un programa en CUDA C es prácticamente código C con algunas funciones especiales referidas a la GPU. A la hora de ejecutar un kernel (lo cual vamos a entender como una función que se ejecuta en la GPU) se debe especificar una serie de parámetros: el número de bloques y el número de hilos. Pero, ¿qué son estos bloques e hilos? Los hilos son la unidad indivisible de CUDA, significando que cada uno se ejecutará de manera independiente en un núcleo de la GPU en un determinado instante. Por otro lado, los bloques no son más que una manera de organizar los hilos tanto a nivel hardware como software. Esta estructura nos permite asignar índices en 3 dimensiones tanto a los hilos como a los bloques. De esta manera, podemos identificar completamente a cualquier hilo perteneciente al conjunto, lo cual facilita en gran medida muchas tareas de programación. Esto se puede ver claramente en la figura 2.2 donde existe una malla bidimensional de bloques y cada uno de ellos tiene a su vez en su interior una matriz de hilos, siendo estos últimos la unidad más pequeña. 6Capítulo 2. Unidades de procesamiento gráfico y CUDA Figura 2.2 Representación esquemática de la jerarquía bloque-hilo 2. Este concepto no se reduce únicamente a 2D sino que es extrapolable a 3D, como puede verse en la figura 2.3. Figura 2.3 Representación tridimensional de una malla 3. Es decir, nosotros somos los que indicamos la cantidad de hilos que se ejecutarán y también la manera en que se organizarán. Esto cobrará mayor importancia cuando expliquemos la diferencia entre programación paralela y estándar (a nivel de código). Cabe aclarar que no se puede generar una cantidad arbitraria de bloques e hilos, sino que esto viene limitado por el modelo concreto de la GPU. Estas limitaciones incluyen un número total de hilos tanto en cantidad total como a lo largo de cada dimensión (que no es estrictamente lo mismo). Otra parte importante de la arquitectura viene representada por los tipos de memoria que existen en su interior. Podemos diferenciar 4 tipos [2]: •La local memory, propia de cada hilo e innaccesible para los demás. 2Imagen procedente de: http://docs.nvidia.com/cuda/cuda-c-programming-guide/#thread-hierarchy 3Imagen procedente de: http://www.karimson.com/posts/introduction-to-cuda/ 2.2 Arquitectura CUDA 7 • La shared memory, la cual es compartida por todos los hilos de un mismo bloque pero no es accesible a los hilos de otro bloque. La ventaja es que resulta tener un acceso más rápido que la memoria local. Este es un buen momento para explicar los mecanismos de sincronización entre hilos. Simplemente, únicamente es posible coordinar a los hilos (deben ser todos ellos) dentro de un mismo bloque. No se permite, por tanto, coordinar bloques ni todos los hilos de todos los bloques. Está limitado a la sincronización dentro de un bloque. • La texture memory, la cual es de sólo lectura y se encuentra optimizada para accesos siguiendo patrones bidimensionales. Es decir, resulta mucho más favorable la consulta de vectores o submatrices que la de un elemento aislado. •La global memory, visible para los hilos de todos los bloques y con un acceso más lento. Una vez hecho este inciso, merece la pena explicar brevemente las diferencias de un programa CUDA C con uno corriente. 3 Ejemplos de programación CUDA C Una vez esbozados cualitativamente los aspectos más generales de la GPU y sus propiedades, se presentan una serie de ejemplos con el objetivo de presentar algunos factores interesantes de la programación en CUDA. 3.1 Cálculo paralelo En este primer ejemplo vamos a tratar de desarrollar las principales diferencias que surgen para un programa que sume 2 vectores basado en: •La programación habitual en C, para ejecutar en la CPU. •La programación en CUDA C, para ejecutar en la GPU. La idea es tratar de esclarecer las diferencias que suponen a nivel de código una u otra opción. Para ello, se va a plantear un programa que suma dos vectores perfectamente conocidos. A modo de resumen, va a a haber principalmente 2 cambios: •Va a ser necesario cambiar la sintaxis de las ”funciones” y la manera en la que estas se invocan. • Es imprescindible plantear los problemas desde otra perspectiva. Muchos de los algoritmos que planteamos en nuestra imaginación los pensamos inherentemente secuenciales. Necesitamos, por tanto, buscar otra manera de resolver los problemas. Figura 3.1 CPU vs GPU 1. 1Imagen procedente de: http://www.e2matrix.com/blog/wp-content/uploads/2014/01/cpu_vs_gpu-1.png 9 16 Capítulo 3. Ejemplos de programación CUDA C function (ya explicaremos más adelante de qué se trata) para distintas cantidades de datos. La prueba ha sido realizada con un procesador Intel Core i7 5820K 3.3GHz y una Nvidia GeForce GTX TITAN Z. El código que se ha ejecutado se comentará posteriormente cuando expliquemos el funcionamiento de un MEX. En el detalle podemos ver como al principio el tiempo de computación es ligeramente inferior en la CPU. Además, para el caso de la GPU, este se mantiene constante. Esto es completamente congruente con lo que hemos explicado con anterioridad. Es decir, ese tiempo se corresponde casi en su totalidad al traspaso de información entre CPU y GPU. Llegados a un determinado instante, este deja de ser representativo. Avanzando en la gráfica, vemos que el tiempo en la CPU empieza a crecer de una manera desorbitada mientras que el de la GPU se mantiene en unos valores realmente bajos. Este tipo de problemas son los que se benefician de usar una GPU, aquellos en los trabajamos con cantidades enormes de información. En caso contrario (cantidades pequeñas o medias), la mayor dificultad de implementación y otras limitaciones provocan que se trate de una opción enrevesada con la que no estamos seguros de obtener beneficios, llegando en algunos casos a empeorar el comportamiento. 4 Matlab Parallel Computing Toolbox Para el desarrollo de este trabajo de fin de grado se ha hecho uso en gran medida del Matlab Parallel Computing Toolbox. Este incluye una serie funciones que permiten trabajar con la GPU directamente en Matlab además de permitir ejecutar programas en C, facilitando en gran medida el movimiento de la información. El Toolbox presenta 2 posibles tipos de paralelización: •Nivel 1: Paralelismo con una única CPU/GPU. •Nivel 2: Paralelismo para varias CPUs/GPUs en una misma máquina. Existe en realidad un tercer nivel basado en el Matlab Distributed Computing Server, pero está orientado a aplicaciones de supercomputación, con gran cantidad de ordenadores trabajando en un mismo problema. Una vez dicho esto, podemos pasar a comentar brevemente algunas de las funcionalidades presentes. 4.1 GpuArray y Gather Esta llamada permite la creación de una variable (incluyendo vectores y matrices) en la memoria de la GPU, para posteriormente trabajar con ellas. Es imprescindible para poder usar algunas de las funciones que veremos más adelante. Gather supone el paso inverso al anterior. Es decir, se encarga de mover las variables de la GPU al workspace de la CPU. 4.2 Bsxfun Esta función de Matlab permite realizar una serie de operaciones binarias elemento a elemento usando como datos dos vectores de entrada. Para ello, es necesario que sea una función del tipo: c=f(a,b) Siendo aybvectores de igual longitud. Los usos que puede recibir esta llamada están relacionados habitualmente con la evaluación de una función para gran cantidad de parámetros distintos. Además, se permite que el propio usuario programe su propia función. Existen varias deficiencias que hacen que no sea de mucha utilidad: •Sólo acepta como parámetros de entrada 2 escalares. Ni siquiera puede prescindirse de uno de ellos. •Únicamente devuelve un escalar. No es posible ampliar la cantidad de variables de salida. •La función no puede acceder a un elemento arbitrario del vector, solo al que le ”corresponde”. Para mostrar su uso, vamos a realizar la resta de dos vectores iguales de 10 componentes. 17 18 Capítulo 4. Matlab Parallel Computing Toolbox Si necesitáramos tener el resultado de c en la CPU no hay más que invocar la llamada gather de la manera: 4.3 Pagefun En el caso de pagefun, se permite aplicar funciones a cada página de una matriz. Matlab define como páginas los valores de la tercera coordenada de una matriz tridimensional. Es decir, estaríamos trabajando con cada una de las matrices bidimensionales que componen la matriz global (figura 4.1). Figura 4.1 Representación gráfica de una matriz tridimensional 1. En este caso, tenemos control total sobre cada ”sub-matriz” , permitiéndonos acceder a cualquier elemento de esta. Además, no existen restricciones en los parámetros de entrada ni de salida por lo que parece una solución más flexible que bsxfun. Sin embargo, tiene una desventaja que provoca que tampoco sea de gran interés y es el hecho de que acepta únicamente unas pocas funciones de Matlab. Esto quiere decir que el usuario no tiene la oportunidad de 1Imagen procedente de: https://www.mathworks.com/help/matlab/math/ch_data_struct4.gif 4.4 Arrayfun 19 programar funciones arbitrarias para usarlas junto a esta llamada. Vamos a ampliar el ejemplo anterior para que muestre la resta de una nube de puntos con otro en un espacio tridimensional. De esta manera, obtenemos en c la resta de los puntos de a (organizados por páginas como vectores columna) con b. Cabe aclarar que el resultado obtenido en ctambién está ordenado por páginas. 4.4 Arrayfun Arrayfun supone una ampliación de lo que era bsxfun. Se adquiere flexibilidad a la hora de decidir los parámetros de entrada y salida. Sin embargo, sigue poseyendo algunas restricciones que tenía su hermana mayor: sólo puede devolver escalares, tomar operandos escalares y no tiene la posibilidad de acceder a otros elementos de los vectores de entrada. Para mostrar el comportamiento de arrayfun se ha programado una función tal que: De manera que la usaríamos así: 20 Capítulo 4. Matlab Parallel Computing Toolbox Podemos ver como hemos ganado en versatilidad con respecto a bsxfun, puesto que nos ha permitido tener 3 variables de entrada y devolver 2 a la salida. Seguimos teniendo aún así la limitación que nos impide devolver otra cosa que no sea un escalar en cada iteración. 4.5 Ejecución de Kernels en CUDA C La ejecución de Kernels desde Matlab es una de las opciones más interesantes que tenemos actualmente debido principalmente a las siguientes razones: • Tenemos una gran flexibilidad a la hora de definir las variables de entrada, salida y las operaciones que se realizan en su interior. Es decir, no estamos restringidos a devolver escalares, aceptar una serie de componentes, etc. Además, tenemos control total sobre la malla y la cantidad de hilos que ejecutan el Kernel. • Es posible decidir cómo gestionar algunos tipos de memoria de la GPU a los que Matlab no era transparente, como la constant memory. •Libera de la complejidad que supone la programación de un fichero MEX (mucho más engorroso). • Debido a que se comportan como funciones, son fácilmente implementables en un fichero MEX que se construya posteriormente. Sin embargo, presentan una serie de inconvenientes: • No permiten la reserva de memoria dinámica durante su ejecución. Por esta razón, es necesaria reservarla anteriormente en Matlab. • A pesar de que los Kernels son realmente una herramienta muy flexible, no somos capaces de solucionar un problema entero con la ejecución de una única función. Esto conlleva a que necesitemos de varios de ellos, con el overhead correspondiente asociado en cada ejecución. Ya hemos visto anteriormente ejemplos de Kernels en C cuando explicamos las diferencias entre la programación en la CPU y la GPU. La diferencia es que ahora lo ejecutamos desde Matlab, por lo que hacemos explícito cómo sería necesario llamarlo desde este entorno: – Es necesario ”inicializarlo”. Es decir, tenemos que definir que archivo lo contiene, la malla y la cantidad de hilos por bloque con la que queremos ejecutarlo. – Una vez hecho esto, se ejecuta a través de la llamada feval. En el primer argumento indicamos el ”objeto” y, posteriormente, van todos las variables de entrada y salida. 4.6 Ejecución de ficheros MEX en CUDA C 21 4.6 Ejecución de ficheros MEX en CUDA C La resolución de problemas mediante la ejecución de ficheros MEX desde Matlab es la mejor solución actualmente. Tenemos una serie de ventajas tales como: •Control total sobre el uso de los diferentes tipos de memoria, variables, etc. •Nos permite acceder a librerías externas como por ejemplo CUBLAS. •Podemos reservar memoria de manera dinámica. El principal inconveniente del uso de este tipo de ficheros es el fuerte incremento de la complejidad con respecto a las anteriores posibilidades. Esto es así en primer lugar por el hecho de que ahora es necesario programar una aplicación entera en lenguaje C, por lo que perdemos las facilidades que Matlab nos proporciona en algunos aspectos. Además, hay que seguir estrictamente una serie de reglas para que todo funcione correctamente. Vamos por tanto a explicar la estructura y funcionamiento de un MEX a partir de uno ya programado. Este ejemplo recibe dos parámetros: Una nube de puntos y la localización de otro punto. El objetivo es buscar en la nube de puntos aquel que se encuentra más cerca de este último. Por tanto, el programa va a tener que realizar una serie de operaciones básicas como: •Definición de los Kernels y librerías. •Inicialización. •Cálculo de la distancia. •Búsqueda del valor mínimo. •Devolución de los resultados. El aumento de complejidad a la hora de programar el MEX se encuentra en la segunda y última etapa, las cuáles están muy relacionadas entre sí. 4.6.1 Definición de los Kernels y librerías Este apartado corresponde exactamente a lo mismo que describimos anteriormente. Los Kernels que invocamos desde Matlab tienen la misma estructura que los correspondientes a los ficheros MEX, únicamente cambia la manera de llamarlos. Aún así, comentamos brevemente la función de cada uno. • La primera función se trata de una invocación en la CPU para calcular el logaritmo en base 2 a partir de las propiedades de los logaritmos. Lo usaremos para decidir la cantidad de bloques que lanzar. • Copia _ vector hace exactamente lo que su nombre indica. En realidad, también podrían usarse llamadas del tipo CudaMemcpy. • El objetivo de Compara es encontrar el índice del punto mínimo una vez este ha sido calculado. Esta operación no se puede realizar directamente debido al algoritmo que se ha usado para el cálculo del mínimo. •La función de la llamada Norma es calcular la raíz cuadrada de la resta. • En Min es donde se realiza el cálculo del mínimo. Además, nos va a dar pie para poner énfasis otra vez en la necesidad de desarrollar nuevos algoritmos. • La llamada Resta realiza la resta de la nube con respecto al punto concreto. Es un paso previo al cálculo de la norma. 24 Capítulo 4. Matlab Parallel Computing Toolbox 4.6.2 Inicialización En la página anterior puede verse el código correspondiente a la inicialización. Aprovechamos para detenernos a revisar la estructura que tiene el inicio de un fichero MEX: • La definición de la mexFunction es siempre la misma. No puede cambiarse el tipo de dato que devuelve, su nombre o los parámetros que recibe. –nlhs. Es un número entero que indica la cantidad de variables de salida. –plhs. Es una tabla que almacena todas las salidas. –nrhs. Es un número entero que indica la cantidad de variables de entrada. –prhs. Es una tabla que almacena todas las entradas. • Definimos una serie de punteros como mxGPUArray. Los que llevan añadido el ”const” corresponderán a valores de entrada y, por tanto, serán constantes. Sin embargo, estos sirven únicamente para comunicarse con Matlab, por lo que necesitamos un puntero accesible desde nuestro programa. El criterio que se ha seguido para los nuevos punteros ha sido añadir d_ a los nombres anteriores. Además, generamos algunas variables auxiliares que usaremos a lo largo del programa. • Una vez hecho esto, es necesario usar la llamada mxInitGPU puesto que es imprescindible para cargar la librería GPU de Matlab. • Una vez finalizado lo anterior, generamos algunas variables dim3 que nos ayudarán a distintos propósitos. Algunas de ellas nos permitirán crear la malla deseada de bloques e hilos. Otras, en cambio, las usaremos para indicar la dimensión de los vectores de salida. • Más adelante, creamos los objetos mxGPUArray y obtenemos los punteros a los datos, tanto de entrada como de salida. En algunas de estas sentencias se hace uso de algunas macros incluidas en las librerías, con el objetivo de inicializar de una manera determinada. •Una vez hemos terminado de declarar variables, pasamos a poner valores en algunas de ellas. •A partir de aquí, ya podemos usar llamadas del tipo cudaMalloc con normalidad. 4.6.3 Cálculo de la distancia Para el cálculo de la distancia haremos uso de varios de los Kernels vistos anteriormente. 4.6.4 Búsqueda del valor mínimo El fragmento de código posterior sintetiza el cálculo del mínimo en la GPU. En el subapartado siguiente nos detenemos brevemente en el método utilizado. 4.6 Ejecución de ficheros MEX en CUDA C 25 Cálculo del Mínimo Para calcular el valor mínimo de un vector en la CPU existen muchas formas de hacerlo. La mayoría de los métodos están basados en recorrer el vector entero y reordenarlo o bien comparar cada elemento con el valor mínimo actual, actualizándolo cuando encontramos uno menor. Todas estas soluciones están basadas en una programación enteramente concurrente, por lo que no es posible trasladarlos a la GPU. La solución propuesta, por tanto, trata de aprovechar las virtudes de la GPU: la gran cantidad de núcleos que posee. El objetivo es lanzar un hilo por cada 2 elementos del vector y comparar cuál de los 2 es inferior, de manera que nos quedaremos únicamente con ese. Este proceso se repite de manera iterativa hasta que sólo quede uno de ellos. Este esquema queda descrito en la figura 4.2. Figura 4.2 Representación gráfica del cálculo del mnínimo.. 32 Capítulo 5. Modelado del sistema cuando esto ocurra, siendo Sla envoltura antes mencionada y ˜xk= [xk,uk]. Para que nos vayamos haciendo una idea, la envoltura convexa implicará que se cumplen una serie de restricciones de igualdad y desigualdad las cuáles expondremos más adelante. A esta serie de puntos que forman una envoltura convexa S (calculados en el paso anterior) los denominaremos estados candidatos xc i∀i=1,...,nc siendo nc el número de estados candidatos, que estará limitado por el número de trayectorias de la base de datos. Ahora entramos en un proceso iterativo en el cual buscamos obtener la envoltura S con el menor número posible de candidatos (ver figura 6.14), con el objetivo de no aumentar la carga computacional en el problema de optimización. (a) ˜xk∈S.(b) ˜xk/∈S. Figura 5.4 Simplificación en R2de la búsqueda de envoltura convexa [12]. De esta manera, variando el valor de ε , aumentaremos o disminuiremos la cantidad de estados candidatos, puesto que cambiarán los puntos que quedan dentro o fuera de la esfera. El método elegido para incrementar o decrementar el valor de ε es el método de la bisectriz, definiéndose: •Partiendo de la resolución del problema para un ε0, se nos plantean dos posibilidades: – Obtenemos una solución factible. Esto significa que puede existir una solución mejor con menor cantidad de puntos, por lo que desplazamos εde la manera que se muestra en la figura 5.5. Figura 5.5 Ejemplo de evolución de epsilon. Los puntos verdes indican que se encontró envoltura convexa. Los puntos azules indican el estado actual, por tanto aún desconocido. 5.2 Modelo basado en datos 33 – No se encontró envoltura convexa. Esto implica que tengo que aumentar el valor de ε para poder incluir más puntos en el rango. Para ello, duplicamos su valor. •Si estoy en una iteración distinta de la primera, existen 4 posibilidades: – Si he encontrado una envoltura convexa y en la iteración anterior se dio el mismo caso, seguiré resolviendo el problema hasta que εk−1−εk≤Tol, tal como se veía en la figura 5.5. – Si no he sido capaz de formar un conjunto S que contenga a ˜xk y sí lo conseguí en el instante anterior, significa que he reducido demasiado la cantidad de puntos de la región. Por tanto, para volver a conseguir una envoltura convexa, necesito avanzar en la dirección contraria tal como se ve en la figura 5.6. Figura 5.6 Ejemplo de evolución de epsilon. Los puntos rojos indican que no se encontró envoltura convexa, al contrario de los verdes. Los puntos azules indican el estado actual, por tanto aún desconocido. – Si anteriormente no fuimos capaces de formar una envoltura convexa y en este instante tampoco, seguiremos duplicando el valor de ε hasta que esto suceda. Esto se ve limitado por otra especificación del algoritmo, εmax , la cual especifica el radio máximo de la esfera de R6 que contiene a los estados. Todo esto se ve reflejado en la figura 5.7. Figura 5.7 Ejemplo de evolución de epsilon. El valor de εmax supone una barrera para el aumento de ε. –Si es la primera vez que encuentro una región convexa, volvemos al primer caso. Sin embargo, seguimos teniendo algunos problemas: • Existen puntos para los que no es posible encontrar una envoltura convexa. Esto puede darse cuando trabajamos con bases de datos demasiado pequeñas o estamos en puntos muy cercanos a los límites del espacio de trabajo del sistema. Para estos casos especiales, se propondrá una alternativa distinta a las restricciones del problema de optimización. • Como ya hemos dicho anteriormente, no nos interesa tener una gran cantidad de candidatos puesto que complican innecesariamente el problema de optimización. Por tanto, decidimos establecer una cantidad de puntos máximo, lo cual rompe un poco la estructura del algoritmo anterior. Sin embargo, nos basta saber que nos proporcionará un mejor rendimiento y un mejor resultado. Esto tiene sentido debido a que una cantidad excesiva de estados candidatos es muy probable que incluya puntos que se encuentran relativamente lejos del estado actual, lo que significa que probablemente no aporten información satisfactoria. Además, aumentan el tamaño del problema a resolver, lo que implica una gran ralentización del proceso de predicción. Resolución del problema de optimización El objetivo de obtener un conjunto S que formara una envoltura convexa (que aún no hemos definido) es poder resolver un problema de minimización con el que obtener unos pesos (a partir de ahora α ) que nos 34 Capítulo 5. Modelado del sistema permitirán calcular una predicción del estado siguiente. Siendo ˜xc+ i el estado siguiente al correspondiente ˜xc i , buscamos estimar el estado tal como se ve reflejado en las ecuaciones 5.9. ˜xk+1= nc ∑ i=1 αi·˜xc+ i(5.9) En función de los resultados obtenidos en apartados anteriores, seguimos un camino u otro: • Hemos conseguido salir del bucle encontrando una envoltura convexa. En ese caso resolveremos el problema con las restricciones genuinas, tal como se ve en las ecuaciones 5.10. Minimizar nc ∑ i=1 α2 i Sujeto a nc ∑ i=1 αi=1(5.10) nc ∑ i=1 αi˜xc i=˜xk∀i=1,...,nc αi≥0∀i=1,...,nc • Si por el contrario, no fuimos capaces de encontrar una envoltura convexa para ε=εmax o bien para N<Nmax , tendremos que liberar alguna restricción para que sea posible obtener alguna solución. En concreto, permitimos que los αi puedan tener valores ≤0 , quedando el problema definido como se refleja en la ecuación 5.11. Minimizar nc ∑ i=1 α2 i Sujeto a nc ∑ i=1 αi=1(5.11) nc ∑ i=1 αi˜xc i=˜xk∀i=1,...,nc −1≤αi≤1∀i=1,...,nc Una vez calculados los pesos, es posible predecir el estado posterior usando la ecuación 5.9. Sin embargo, antes de continuar, vamos a detenernos en diversas posibilidades para la resolución de este problema que acabamos de plantear. Programación cuadrática Existen multitud de algoritmos dedicados a resolver problemas de optimización. Dependiendo del software que usemos, tendremos unas posibilidades u otras. En nuestro caso, al estar tratando de realizar la mayor cantidad posible de implementaciones en la GPU, nos planteamos la posibilidad de resolver también este problema en su interior, además de los anteriormente expuestos. Por ello, nos basamos en la literatura que hay actualmente sobre algoritmos de programación cuadrática que fueran susceptibles de ser paralelizables. De esta manera, encontramos información en diversos artículos como ”A Parallel Quadratic Programming Algorithm for Model Predictive Control” [13], ”Parallel Quadratic Programming for Image Processing” [14] y ”Projection-free Parallel Quadratic Programming for Linear 5.2 Modelo basado en datos 35 Model Predictive Control” [15]. En ellos se desarrolla un algoritmo completamente paralelizable que se puede programar en la GPU. En algunos de los trabajos anteriores se usa para resolver los problemas que aparecen en las técnicas de control predictivo basado en modelo, mientras que en otros, tiene como finalidad la mejora computacional de programas de tratamiento de imágenes. Partiendo de un problema de optimización como el 5.12, tenemos unas variables de decisión Uk con una matriz Q que multiplica a los términos cuadráticos y otra H para los lineales. Tenemos a su vez una serie de restricciones de desigualdad representadas por las matrices VyW. Minimizar 1 2UTQU +HTU Sujeto a (5.12) VU <W Sin embargo, el algoritmo no es capaz de trabajar directamente con problemas expresados de esta manera. Por tanto, es necesario transformar las ecuaciones anteriores a su forma dual, obteniendo como resultado el problema 5.13. Minimizar 1 2yTqy +hTy Sujeto a (5.13) y>0 Las nuevas matrices qyhquedan definidas por las ecuaciones 5.14 y 5.15. q=VQ−1VT(5.14) h=W+VQ−1H(5.15) Una vez resuelto el problema, es necesario calcular la solución U∗ a partir del óptimo y∗ siguiendo la expresión 5.16. U∗=−Q−1(H+VTy∗)(5.16) Sin embargo, todavía no podemos aplicar el algoritmo a nuestro problema concreto. Esto es debido a que no hemos tenido en cuenta las restricciones de igualdad. Se plantearon en un principio diversas posibilidades como el uso de ”penalty functions”. Finalmente, se decidió despejar un conjunto de variables igual al número de restricciones, de manera que se obtiene un sistema compatible indeterminado. Para resolver este sistema, se hace uso de la factorización QR [16]. Antes de presentar los pasos que se han de dar, conviene mostrar qué representa cada variable. A=Matriz que aglutina las restricciones de desigualdad b=Vector que contiene el valor de las desigualdades Q,R=Factorización QR de la matriz AT p=Número de restricciones a eliminar Las otras matrices Q1 y Q2 provienen de la misma matriz Q . Esta subdivisión es dependiente del número de restricciones de igualdad que queramos transformar. Dicho esto, se procede de la siguiente manera: 36 Capítulo 5. Modelado del sistema •Se calcula la factorización QR. Eso implica que hemos obtenido QyR. •Dividimos las matrices de la siguiente manera: QT 1Q1=ISiendo Q1una matriz que contiene las primeras pcolumnas de la matriz Q. QT 2Q2=ISiendo Q2una matriz que contiene el resto de columnas de la matriz Q. Se cumple asimismo que QT 1Q2=0 . Esto puede servir para comprobar que las operaciones se han realizado correctamente. •Una vez hecho esto, es posible calcular una de las infinitas soluciones del sistema con la expresión: ˆx=Q1R−Tb • LLegados a este punto, es conveniente recordar el problema 5.12. Vamos a nombrar a partir de ahora Qeq y Heq a las matrices que resultan de la eliminación de las restricciones de igualdad. Es posible calcularlas a partir de las expresiones siguientes: Qeq =QT 2QQ2 heq = (hT+ˆxTQ)Q2 • Además, es necesario modificar las restricciones de desigualdad, puesto que también se ven afectadas por la transformación que acabamos de realizar. Por tanto, siguiendo el mismo criterio anterior, definiremos Veq yWeq de la siguiente manera. Veq =−Q2 Weq =ˆx+W • Por último, una vez se hubiera resuelto el problema, es necesario revertir la transformación (de manera parecida a lo que ocurría con el problema dual). Para ello, simplemente hacemos uso de: U∗=Q2x∗+ˆx Sin embargo, nuestro problema original (5.10) posee tanto restricciones de igualdad como de desigualdad. Por tanto, necesitamos aplicar en primer lugar la eliminación de las primeras y, una vez hecho esto, convertirlo al problema dual. Esto va a requerir de igual manera cuando obtengamos la solución que haya que revertir estas transformaciones. Conviene tener en mente todas las operaciones que hemos tenido que realizar para llegar a este punto, puesto que suponen un punto clave a la hora de la resolución que hemos tomado. El algoritmo que pasamos a desarrollar es llamado Parallel Quadratic Programming y es comúnmente conocido por sus siglas (PQP). Sus iteraciones están basadas en la realización de operaciones matriciales sencillas las cuáles pueden ser fácilmente paralelizables. Esto no quiere decir que únicamente se obtenga beneficio en el caso de la implementación en GPU, sino que es capaz de equipararse a muy diversos algoritmos ysolvers comerciales [13]. Los pasos a seguir por el algoritmo son los siguientes: 1. A partir de la formulación genérica de programación cuadrática, eliminar las restricciones de igualdad usando, por ejemplo, la factorización QR. 2. Una vez definido este nuevo problema derivado del anterior, reescribirlo a su forma dual correspondiente. 5.2 Modelo basado en datos 37 3. Empezando con un valor inicial y0>0 , aplicar la siguiente actualización hasta alcanzar la tolerancia especificada. yi←yih− i+ (Q−y)i h+ i+ (Q+y)i,(5.17) Siendo: q+=max(q,0) + diag(r) q−=max(−q,0) + diag(r) h+=max(h,0) h−=max(−h,0) No se ofrece mucha información en los artículos anteriores respecto al vector r . Nos dicen que han comprobado que funciona para valores cercanos a 0, pero no especifican cuál es exactamente su función o cómo afecta al algoritmo cambiar su valor. Si nos fijamos bien en la ecuación 5.17, podemos ver que dentro de los corchetes se obtendría tanto arriba como abajo un vector. En realidad, esta división se realiza elemento a elemento, posteriormente actualizando a la variable de decisión correspondiente. Por esta razón es posible paralelizar el algoritmo, porque la actualización de cada una de las variables de decisión se realiza por separado. Se planteó más adelante una posible mejora para el rendimiento del algoritmo invocando una actualización distinta cada cierto número de ejecuciones [15]. Sin embargo, no resultaba ser muy estable en algunas circunstancias. En algunos casos conseguía mejorar los tiempos mientras que en otros los empeoraba. 4. Revertir la solución a su forma primal. Una vez hecho esto, devolverlo a su forma primera. Es importante seguir el orden aquí especificado. Esto es debido a que hay que ser congruente con el método que se siguió a la hora de transformarlos en primer lugar. Es fácil darse cuenta de que se sigue el orden contrario al que se siguió a la hora de la conversión (figura 5.8). Figura 5.8 Conversión de restricciones del problema de optimización. 38 Capítulo 5. Modelado del sistema Como podemos ver, la carga computacional del algoritmo PQP es realmente pequeña. Sin embargo, nuestro problema concreto necesita pasar por 2 transformaciones consecutivas antes de poder resolverse. Además, algunas de las operaciones involucradas conllevan inversiones de matrices, las cuáles suelen ser computacionalmente costosas y en nuestro caso tendrían que realizarse para cada iteración del método de la bisectriz. Por tanto, al tratar de comparar su rendimiento con el solver de Matlab (quadprog), hemos preferido este último por diversas razones: • Quadprog resulta ser mucho más maduro que el algoritmo anteriormente presentado y, por tanto, no diverge para algunos casos patológicos como ocurría con el PQP. • A pesar de que el PQP era generalmente más rápido que quadprog, las operaciones necesarias para devolver los resultados incrementaban este tiempo por encima del obtenido directamente con el solver de Matlab. • Aunque quizá no parezca una razón al uso, el código se vuelve mucho más limpio y no es necesario generar una serie de variables auxiliares que estarían consumiendo memoria en la GPU. Aún así, existen problemas en los que nos podríamos beneficiar de las limitaciones del PQP, como aquellos en los que la conversión del problema dual pueda realizarse fuera de línea debido a que permanece invariable. En esos casos, no tendríamos la carga computacional asociada a las inversiones de matrices sino únicamente la ejecución del algoritmo que, como hemos dicho antes, suele ser más rápida que quadprog. 5.2.3 Resultados De igual manera que la bondad de los resultados que se obtienen al modelar un determinado sistema a través de una función de transferencia depende de factores muy diversos como el orden del modelo, las dinámicas que no se tienen en cuenta por simplicidad, etc, existen en nuestro caso factores que favorecen o empeoran el comportamiento. Algunos de estos parámetros serán: • En primer lugar, como es obvio, la cantidad de información que tenemos en la base de datos. Es evidente la razón por la que esto es así, no es posible predecir el estado de un determinado sistema si no conozco mínimamente como ha evolucionado anteriormente. Además de la cantidad, es importante la selección que hemos hecho, es decir, qué información hemos decidido almacenar. Se puede dar el caso de tener una base de datos enorme, pero que no está muestreada correctamente (trayectorias repetidas o redundantes o mala selección de las trayectorias a almacenar en general), de manera que no es capaz de predecir adecuadamente. • El εmax . Si no existiera Nmax , εmax afectaría muchísimo al rendimiento al aumentar su valor, puesto que no habría límite de estados candidatos. Diciéndolo de esta manera parece que no tendría sentido tener un valor grande de εmax , esto no siempre es así. En función de la zona en la que estemos trabajando y de como esté construida la base de datos, podemos encontrarnos para una misma trayectoria con que para un mismo valor de εmax y distintos instantes de tiempo, tenemos algunos casos en los que no se encuentra ningún candidato (posible en los límites del espacio de trabajo) mientras que en otros se supera el máximo (esto ocurre sobre todo cuando nos acercamos al régimen permanente). •Nmax . Se encarga principalmente de impedir que se acumulen demasiados estados candidatos, de manera que impide que se pueda perder demasiado tiempo a la hora de resolver el problema de programación cuadrática. Se complementa bien con εmax , puesto que ahora nos podemos permitir aumentar su valor sin miedo a que la ejecución se haga muy lenta. Para comprobar el funcionamiento del modelo, se han generado trayectorias que parten de estados iniciales aleatorios con referencias igualmente aleatorias. En las figuras 5.9 y 5.12 puede verse el funcionamiento del modelo con una base de datos suficientemente amplia. El error de la predicción es realmente pequeño, por lo que se puede combinar con el uso de sensores industriales sin problema. El funcionamiento para el caso de la figura 5.10 resulta hasta sorprendente. Hemos reducido la cantidad de trayectorias hasta 1000, es decir, 100 veces menos que en el anterior. Aún así, seguimos obteniendo resultados 5.2 Modelo basado en datos 39 bastante satisfactorios. Se ha realizado además la predicción de un conjunto de 1000 trayectorias para los dos conjuntos de bases de datos. La comparativa de errores queda reflejada en la tabla 5.2. Tabla 5.2 Comparativa de la bondad para diferentes tamaños de la base de datos. Nºtrayectorias Error máximo Error mínimo Error medio 1000 0.0336 2.8228 ·10−60.0019 100000 0.0195 0 4.336 ·10−4 Como era de esperar, los resultados son mejores para una base de datos más grande. Sin embargo, para la precisión que tienen gran cantidad de sensores industriales, el conjunto de 1000 trayectorias funcionaría sin problemas. 40 Capítulo 5. Modelado del sistema Figura 5.9 Modelo con una base de datos que almacena 100000 trayectorias. Ejemplo sin ruido. 5.2 Modelo basado en datos 41 Figura 5.10 Modelo con una base de datos que almacena 1000 trayectorias. Ejemplo sin ruido. 48 Capítulo 6. Control predictivo Lo que planteamos a partir de ahora se complementa perfectamente con el modelo visto anteriormente, cuyo algoritmo presenta ciertas similitudes. Partiendo del siguiente modelo no lineal discreto: xk+1=f(xk,uk)(6.3) Donde xk∈Rn es la salida del sistema y uk∈Rm es la entrada en el tiempo de muestreo k respectivamente. frepresenta el modelo del sistema, el cual se supone desconocido. Queremos diseñar una ley de control que permita resolver el problema de regulación de manera que podamos partir de un x0arbitrario para alcanzar una referencia xre f cualquiera (figura 6.4). Figura 6.4 Problema de regulación en R2. Buscamos hallar la entrada actual ukde la manera: uk=DbPC(xk,xre f ,BD,ID)(6.4) Los ID (Índices de desempeño) representan la bondad de las trayectorias candidatas y están definidos de la forma siguiente: IDk= N ∑ i=0 Ji(h,q,c,p)(6.5) Siendo IDk el índice de desempeño correspondiente a una trayectoria candidata k , N el horizonte de predicción y Jiel coste de etapa. La definición de Jiviene dada por: Ji(h,q,c,p) = q2 a i +c·q2 b i +pVmin A(h1i+h2i)∀i=1,...,N(6.6) Los valores cypson grados de libertad que nos permiten dar un mayor peso a unos u otros términos. 6.2 Controladores predictivos basados en datos 49 Estos índices de desempeño no se encuentran almacenados en la base de datos, sino que se calculan en línea. Sin embargo, no se trata de una operación computacionalmente dominante, ni siquiera en la CPU. Si nos remontamos al apartado anterior donde explicábamos el MPC, es fácil darse cuenta de la similitud de lo que suponen estos índices con la función de costes cuya optimización nos resultaba de gran importancia para obtener las secuencias de control óptimas. 6.2.1 Algoritmo de control A pesar de ser realmente parecido al usado anteriormente para el modelado, conviene pararse brevemente en algunos aspectos concretos del algoritmo. El esquema viene representado en la figura 6.5. Figura 6.5 Esquema del algoritmo de control basado en datos. • La búsqueda de estados candidatos es prácticamente la misma con ciertas salvedades. Antes se trataba de una búsqueda en R6 con un vector del tipo ˜xk= [xk,uk] . En el caso anterior, se añadían las uk al 50 Capítulo 6. Control predictivo modelo para poder tener información de la tendencia del sistema. Sin embargo, para el control, lo que nos interesa es conocer el punto objetivo, puesto que con el estado actual y las referencias, tenemos toda la información necesaria. Eso quiere decir que tendremos ˜xk= [xk,hre f ]. Por otro lado, antes sólo era necesario tener un estado siguiente para poder predecir. En nuestro caso actual, esto dependerá del horizonte de predicción. Esto quiere decir que ahora la candidatura está más restringida. • Una vez hemos calculados los candidatos, se comprueba que existe una envoltura convexa. Esto no quiere decir que ya estemos calculando una posible solución al problema de optimización sino que simplemente comprobamos la factibilidad. Además, las restricciones que venían impuestas anteriormente han cambiado (junto al problema de optimización). Las nuevas restricciones vienen definidas por las siguientes ecuaciones: nc ∑ i=1 α=1∀i=1,2,...,nc (6.7) nc ∑ i=1 αixc i=xk∀i=1,2,...,nc (6.8) ∑αixre f i=xre f ∀i=1,2,...,nc (6.9) αi≥0∀i=1,2,...,nc (6.10) A la vista de las ecuaciones anteriores, vemos que el problema se descompone en la búsqueda de dos envolturas: una de R4 que envuelva las alturas de los tanques y otra de R2 que cobije a las referencias. Una simplificación gráfica de esto puede verse en la figura 6.6. Figura 6.6 Ejemplo simplificado de búsqueda de envoltura convexa para el problema de control. Caso de solución factible encontrada. Si no se encuentra solución factible al conjunto de ecuaciones anterior (figura 6.7), será necesario repetir la búsqueda anterior con un ε mayor. Además, seguiremos usando el método de la bisectriz 6.2 Controladores predictivos basados en datos 51 para calcular los valores de ε para la siguiente iteración. De igual manera, aplicamos igualmente los criterios de εmax yNmax, con el mismo objetivo que en casos anteriores. Figura 6.7 Ejemplo simplificado de búsqueda de envoltura convexa para el problema de control. Caso de solución factible no encontrada. • El paso siguiente es completamente nuevo y ha sido esbozado previamente. Se trata del cálculo de los índices de desempeño. Estos serán parte imprescindible del problema de optimización que es inherente al control predictivo. •El nuevo problema de minimización viene definido por las ecuaciones siguientes: Minimizar nc ∑ i=1 αiIDi Sujeto a nc ∑ i=1 α=1∀i=1,2,...,nc nc ∑ i=1 αixc i=xk∀i=1,2,...,nc (6.11) ∑αixre f i=xre f ∀i=1,2,...,nc αi≥0∀i=1,2,...,nc Vemos como hemos cambiado el problema anterior por uno de programación lineal donde están incluidos los índices de desempeño que hemos calculado en el paso anterior. De aquí obtendremos los pesos necesarios que nos permitirán calcular la ukóptima para este instante de muestreo. Sin embargo, existe un caso que no hemos estudiado, el correspondiente a que no existe una envoltura convexa para cualquier valor de ε que satisfaga ε≤εmax y N≤Nmax . Para los puntos en los que esto suceda, hemos decidido quedarnos con la acción de control correspondiente al punto más cercano (según la distancia R6previamente definida). 52 Capítulo 6. Control predictivo 6.2.2 Tareas en la GPU Gran parte del objetivo de este proyecto es aplicar la GPU en posibles problemas de control pero, aún no hemos hablado mucho del papel que esta realiza tanto en el modelo como en el control. En los algoritmos que hemos presentado anteriormente diferenciamos básicamente los siguientes problemas: • Búsqueda de candidatos. Este problema incluye el cálculo de la distancia de una nube de puntos a otro y el cálculo de los mínimos de cada trayectoria. • Optimización de una determinada función para hallar unos pesos α . Incluimos también aquí la comprobación de las restricciones. •Cálculo de los índices de desempeño. Resulta claro que la resolución de problemas de minimización no se va a ver recompensada por el uso de la GPU debido a lo que comentábamos del PQP en el capítulo anterior. Además, para que este se viera mejorado, necesitaríamos problemas con gran cantidad de variables de decisión, lo cual no es nuestro caso. Por otro lado, el cálculo de los índices de desempeño es un problema muy fácilmente paralelizable debido a que el ID de cada trayectoria es independiente del de la siguiente. Sin embargo, debido a que la cantidad de IDs a calcular no es suficientemente grande y que las operaciones que se realizan son realmente sencillas, no se experimenta mejoría en este aspecto, pero tampoco empeora. Igualmente, en el código final se encuentra programado en la GPU. Por tanto, es fácil ver que el problema que se va a ver más beneficiado de la paralelización va a ser la búsqueda de los candidatos. Esto va a ser así principalmente por el gran volumen de información que hay que procesar. Es decir, para una base de datos ”pequeña” de 1000 trayectorias, tenemos aproximadamente unas 2 millones de filas. Por tanto, eso supone que hay que calcular la distancia del estado actual a unos 2 millones de puntos. Es fácil darse cuenta de que si esto es así para el caso ”pequeño”, el grande tendrá un volumen mucho mayor, lo que llevará a un mejor rendimiento en comparación a la CPU. 6.2.3 Resultados de simulación Vamos a estudiar distintas características del controlador en función de los tamaños de las bases de datos, como la evolución de los índices de desempeño, la factibilidad del problema de optimización, el coste computacional y la mejora con respecto a la CPU. Para obtener estos resultados, se han realizado ensayos con las referencias propuestas en el concurso de la CEA del año 2014. Además, se ha implementado sobre un modelo realista de la planta, de manera que presenta ruido en la medida. Las figuras siguientes muestran los resultados de ambos ensayos (con ruido y sin él) para el caso de 100000 trayectorias. 6.2 Controladores predictivos basados en datos 53 Figura 6.8 Resultados de la simulación con un modelo sin ruido. 54 Capítulo 6. Control predictivo Figura 6.9 Resultados de la simulación con un modelo ruidoso. 6.2 Controladores predictivos basados en datos 55 Índices de desempeño La evolución de los índices de desempeño puede verse en el gráfico de barras de la figura 6.10. Figura 6.10 Evolución del valor medio de los índices de desempeño. Debido a la sencillez del problema, aumentar el tamaño de la base de datos no mejora en gran medida el rendimiento del controlador, sino que permanece oscilando en torno a un punto. Sin embargo, la recta de regresión de los datos proporciona una pendiente decreciente, lo cual nos podría estar indicando una mejoría. Factibilidad La aparición de puntos no factibles tiene lugar cuando la base de datos no está bien muestreada o bien no tiene información suficiente. Como vemos en la figura 6.11, existe siempre una solución factible al problema de optimización llegados a un determinado tamaño. Figura 6.11 Porcentaje de puntos factibles. 56 Capítulo 6. Control predictivo Coste computacional y mejora Vamos a centrarnos en el estudio del coste de la búsqueda de estados candidatos. Esto incluye tanto el problema del cálculo de la distancia como el de la búsqueda de los mínimos. Además, veremos la mejora respecto a su versión no paralelizada. (a) Mínimo. (b) Distancia. Figura 6.12 Coste computacional (GPU). Como era de esperar, el tiempo necesario para realizar los cálculos aumenta conforme se añaden trayectorias a la base de datos. En realidad, la tendencia debería ser en ambos problemas lineal pero, debido al método de la bisectriz (cuyas iteraciones dependen del número de trayectorias), la cantidad de veces que se ejecuta la búsqueda de los mínimos no es constante, provocando por tanto que el cálculo del mínimo no presente una tendencia lineal. (a) Mínimo. (b) Distancia. Figura 6.13 Mejora de velocidad. Puede verse que la tendencia de estas curvas no es lineal. Esto es debido a que al principio se experimenta una mejora muy grande pero, conforme vamos avanzando esta comienza a mantenerse. Debería llegarse a un punto en el que la GPU sature y que, por tanto, deje de incrementar la mejoría. Aunque no se hayan representado los tiempos de todas las funciones que actúan en el programa, cabe aclarar que la ejecución del programa es lo suficientemente rápida como para cumplir con el tiempo de muestreo, tanto en la versión de la GPU como de la CPU. Es más, las gráficas anteriores se refieren a 1500 instantes de muestreo. Es fácil ver que se supera con creces la imposición temporal de 5 segundos. 6.2 Controladores predictivos basados en datos 57 Hebras y bloques Por otro lado, presentamos cómo afecta la cantidad de hebras por bloque a la hora de la ejecución de los kernels. En un principio, con 1 hilo por bloque, la ejecución se hace realmente lenta y posteriormente cae en picado hasta permanecer relativamente constante, lo cual es completamente congruente con la arquitectura de las GPUs y su funcionamiento. (a) Mínimo. (b) Distancia. Figura 6.14 Tiempo normalizado de computación respecto a la cantidad de hebras por bloque. Cuando nos referimos a tiempo normalizado queremos decir que hemos dividido todos los valores por el máximo de todos ellos, con lo que obtenemos un resultado adimensional. En nuestro caso, hemos elegido una cantidad de 500 hebras por bloque para las pruebas anteriores. 6.2.4 Pruebas en la planta real Una vez visto el funcionamiento en simulación, nos planteamos la posibilidad de implementar el controlador anterior en la planta de los 4 tanques presente en los laboratorios de la escuela. El envío de las señales de control, medidas de los sensores, etc, se hace a partir de la conexión de un ordenador y un autómata programable. Esta comunicación se lleva a cabo a través del protocolo OPC, aunque resulta transparente al usuario. Una vez dicho esto, cabe aclarar que para poder controlar la planta existen 2 posibilidades: •Uso del entorno labview. •Uso de Matlab y Simulink. En nuestro caso, al tener ya definido gran parte del código en Matlab, nos decantamos por este último. Simplemente sería necesario adaptar nuestro fichero ”.m” a un bloque de Simulink. Esto es así porque únicamente se necesita conectar las salidas de lo que sería nuestro controlador ( qa y qb ) a la entrada del bloque OPC encargado de la comunicación, y lo contrario para el caso de las entradas (h). Sin embargo, existe un problema, este no es más que el hecho de que el ordenador que se encarga de realizar estas tareas no disponía de tarjeta gráfica dedicada con arquitectura CUDA. Por ello, usamos un ordenador portátil (que sí dispone de CUDA) el cual conectamos usando un cable ethernet cruzado al ordenador encargado de gestionar la planta, creando una pequeña red de área local (LAN). Pero claro, es necesario transmitir la información de un Matlab a otro, puesto que nuestro objetivo es que el ordenador ”remoto” sea el encargado de calcular las acciones de control usando el algoritmo anterior. Esto quiere decir que en cada instante de muestreo debe enviarse el estado del sistema y sus referencias y, una vez realizados los cálculos, deben devolverse las acciones de control que deben aplicarse al sistema. 64 Índice de Figuras 6.15 Esquema de conexiones de la LAN 58 6.16 Ejemplo de comunicación 59 6.17 Ensayo en la planta real 59 6.18 Simulación con modelo discrepante 60 Índice de Tablas 5.1 Ejemplo de celdas en la base de datos 30 5.2 Comparativa de la bondad para diferentes tamaños de la base de datos 39 65 Bibliografía [1] NVIDIA Corporation. http://www.nvidia.es/object/cuda-parallel-computing-es.html , 14 de Junio de 2017. [2] J. Sanders yE. Kandrot,CUDA by example, 2007. [3] NVIDIA Corporation. http://www.nvidia.com/object/uavs-drones-technology.html , 14 de Junio de 2017. [4] MathWorks. https://es.mathworks.com/products/parallel-computing.html , 14 de Junio de 2017. [5] Grupo de estimación, predicción, optimización y control,Transparencias de programación en CUDA con Matlab. [6] MathWorks. https://es.mathworks.com/help/distcomp/gpuarray.html , 14 de Junio de 2017. [7] MathWorks. https://es.mathworks.com/help/distcomp/gather.html , 14 de Junio de 2017. [8] MathWorks. https://es.mathworks.com/help/matlab/ref/bsxfun.html , 14 de Junio de 2017. [9] MathWorks. https://es.mathworks.com/help/distcomp/pagefun.html , 14 de Junio de 2017. [10] MathWorks. https://es.mathworks.com/help/matlab/ref/arrayfun.html , 14 de Junio de 2017. [11] K. H. Johansson,”The quadruple-tank process”, IEEE Transactions on Control Systems Technology, vol 8 (2000). [12] V. S. Gómez,D. Muñoz de la Peña yT. Álamo,”Modelo basado en optimización y bases de datos”, 2015. [13] M. Brand,V. Shilpiekandula,C. Yao yS.A. Bortoff,”A Parallel Quadratic Programming Algorithm for Model Predictive Control”, Mitsubishi Electric Research Laboratories, Agosto de 2011. [14] M. Brand yD. Chen,”Parallel Quadratic Programming for Image Processing”, Mitsubishi Electric Research Laboratories, Septiembre de 2011. [15] S. Di Cairano,M. Brand yS.A. Bortoff,”Projection-free Parallel Quadratic Programming for Linear Model Predictive Control”, International Journal of Control (2013), 86:8, 1367-1385. [16] S. Boyd yL. Vandenberghe,Convex Optimization, Séptima edición, 2009. [17] D. Rodríguez yC. Bordóns,Apuntes de Ingeniería de Control, Revisión del año 2007. [18] V. S. Gómez,”Control predictivo basado en datos. Aplicación a un sistema de cuatro tanques.” TFG , 2015. [19] V. S. Gómez,D. Muñoz de la Peña yT. Álamo,”Control predictivo basado en datos”, 2015. 67