Colocación óptima de UAVs para la reconstrucción del campo de velocidades del viento en zonas de inspección

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 Tecnologías Industriales Colocación óptima de UAVs para la reconstrucción del campo de velocidades del viento en zonas de inspección Autor: Fernando González Balbontín Tutores: Antonio González Morgado Guillermo Heredia Benot Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Colocación óptima de UAVs para la reconstrucción del campo de velocidades del viento en zonas de inspección Autor: Fernando González Balbontín Tutores: Antonio González Morgado Doctorando Guillermo Heredia Benot Catedrático de Universidad Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Trabajo Fin de Grado: Colocación óptima de UAVs para la reconstrucción del campo de velocidades del viento en zonas de inspección Autor: Fernando González Balbontín Tutores: Antonio González Morgado Guillermo Heredia Benot 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 Ap rovecho para expresar mi agradecimiento hacia todas las personas que han servido de apoyo y me han ayudado, cada una a su manera, en la ejecución de este Trabajo Fin de Grado. Entre todas las personas quisiera destacar en primer lugar a mi tutor, Antonio, que me ha brindado de una excelente guía académica y su apoyo en los momentos más complicados. Siento gran gratitud por haber compartido horas de trabajo con una gran persona, algo que ha impulsado aún más mis ganas de trabajar. Quisiera también expresar mi agradecimiento a todos los profesores que son artífices de lo que soy hoy. Son ellos los que han aportado todos esos conceptos sustanciales que me han permitido llegar hasta aquí. Especialmente me gustaría agradecer a aquellos que han sobrepasado la barrera de los estudios y han sabido tratar temas personales. Por último, me gustaría agradecer a mi familia y mis amigos. Son ellos los que me motivan a seguir creciendo, tanto académicamente como en lo personal, algo por lo que siempre estaré infinitamente agradecido. Fernando González Balbontín Estudiante de Grado en Ingeniería de Tecnologías Industriales Sevilla, 2023 I Resumen Es te estudio lleva a cabo la investigación de la colocación de manera eficiente de sensores instalados en UAVs para llevar a cabo la reconstrucción del campo de velocidades del aire en una zona donde se quiera llevar a cabo la inspección con un UAV. Esta reconstrucción servirá como predicción de viento para el UAV inspector, que podrá modificar la planificación de trayectorias en tiempo real. Para lograr esto, es necesario tener claro ciertos conceptos matemáticos de optimización y estadística, pero también un gran manejo con las ecuaciones de Navier-Stokes para llevar a cabo unas buenas simulaciones. Juntando estas simulaciones con los conceptos matemáticos, se logra efectuar la reconstrucción del campo fluido. La optimización de sensores se basa en la posibilidad de obtener de manera directa datos comprimidos y resolver la alta dimensión de los datos en función de unos pocos representativos. Siguiendo esta línea aparece el concepto sparse. Las aplicaciones en las que puede aplicarse este concepto, son diversas. En primer lugar, aquellas situaciones en las que existen multitud de sensores y se pretende reducirlos, sin dejar de perder información. En segundo lugar, en situaciones en las que los sensores son costosos y se pretende implementar el mínimo número de ellos. En el caso de la reconstrucción del campo fluido, se pretende reducir el coste de la medición lo máximo posible, sin dejar de obtener una información precisa del fluido en el emplazamiento. Una vez determinado el número y la ubicación de los UAVs, se toman las medidas y se realiza la reconstrucción. Esta información valiosa sirve a la hora de tomar decisiones en el terreno de la inspección. III 2Capítulo 1. Introducción algoritmo que hay que usar, es necesario crear un dataset que contenga multitud de simulaciones del campo de velocidades. Se deben plantear tanto los distintos solvers de fluidos para determinar aquel que cumple los requerimientos de calidad, como las diversas simulaciones que han de llevarse a cabo, conteniendo una rica diversidad de distintos ángulos de ataque y velocidades. Por último, una vez que se tenga generado el dataset, se le aplicarán los distintos algoritmos de reducción de dimensionalidad y de reconstrucción de datos, generando unos resultados que determinarán, qué ubicación óptima deben tomar los pequeños UAVs, que se encargan de tomar las medidas de velocidad y orientación de viento, para que el UAV que realiza la inspección genere la predicción de viento. 1.2 Estructura La estructura que ha tomado el proyecto tiene una división en 5 capítulos: • En el Capítulo 1, se realiza una introducción sobre lo importante que es tener en cuenta el viento cuando se hace la planificación de trayectorias de un UAV, y las posibles soluciones. • En el Capítulo 2, se introducen todos los conceptos matemáticos necesarios para el entendimiento del concepto Compressed Sensing. • En el Capítulo 3, se lleva el concepto Compressed Sensing a problemas más focalizados. Con esto, se busca reducir la dimensionalidad del problema para llevar a cabo la búsqueda de sensores ubicados óptimamente, finalizando en una reconstrucción. • En el Capítulo 4, se implementan los algoritmos vistos en capítulos anteriores a un dataset generado en este capítulo. Se obtienen los resultados acerca de donde se tienen que ubicar los UAVs para tomar las mediciones óptimas. • En el Capítulo 5, se analiza brevemente que se ha hecho en todo el trabajo, viendo a su vez las posibles futuras líneas de trabajo. 2 Compressed Sensing En este capítulo se introduce el concepto de Compressed Sensing. El término Compressed Sensing, hace referencia al campo de las matemáticas encargado de tratar la adquisición, tratamiento y reconstrucción de señales a partir de un número reducido de mediciones, menos de las esperadas en un primer momento. Esto último, se basa en el hecho de que la mayoría de señales son sparse cuando son expresados en una base determinada. El concepto sparse, viene a referirse en este contexto que solo unos pocos parámetros son necesarios para caracterizar un conjunto de datos. En concreto, nos informa de qué datos tienen más relevancia y en que proporción, frente otros. Toda compresión de datos lleva de manera inherente este concepto, lo que permite adquirir datos eficientemente, al mismo tiempo que se mantiene la veracidad. Se verán todos los conceptos matemáticos que están detrás de este campo revolucionario, las útiles aplicaciones y los ejemplos directamente relacionados. Como se verá, las áreas de las matemáticas que se trabajarán son el álgebra, el cálculo, la geometría y la probabilidad. Juntando un poco de todas estas disciplinas, es posible llegar a conclusiones bastante potentes como es el Compressed Sensing, tanto en algoritmos eficientes, como en lo que significa llegar a una solución de tal magnitud. 2.1 Sparsity En consonancia con lo comentado anteriormente, la mayoría de señales con las que se trabajan son alta y fácilmente compresibles. La idea que subyace de esta compresibilidad es la siguiente. Se tienen unos datos que están expresados en unos ejes determinados. Lo que se busca, por tanto, es una transformación a partir de la cual, los datos expresados en los nuevos ejes se representen de tal manera que tan solo algunos pocos de ellos tengan relevancia. Esto implica que si se es capaz de guardar solo aquellos modos más enérgicos o, lo que es lo mismo, aquellos datos importantes, podemos eliminar gran cantidad de información que no es relevante. Por tanto, lo que se busca son unos ejes universales, como puede ser la transformada de Fourier, en los que la señal puede ser representada con vector sparse. Centrándonos en terminología matemática, podemos expresar la señal como un vector x∈Rn como otro vector sparse,s∈Rn, generado tras la transformación universal ψ∈Rnxn : x=ψs(2.1) El vector s es el que venimos comentando que tendrá P distintos de cero. La condición para que el vector s sea sparse, se cumple cuando la base ψ es genérica. Por tanto, esos ciertos P términos proyectados en la base genérica ψ , reconstruyen la señal x . Se tendrá que introducir, por tanto, a la base universal por excelencia, la transformada de Fourier [12]. 3 4Capítulo 2. Compressed Sensing 2.2 Fourier y Wavelet Transforms Es posible que la transformada de Fourier sea la más transcendental, versátil y usada transformación hasta día de hoy. Fue J.-B. Joseph Fourier quien introdujo el concepto de transformar señales que están en dominio del tiempo a señales que están en el dominio de la frecuencia. A la conclusión que llegó fue, que una función puede ser expresada como suma de cosenos y senos de distintas frecuencias, formando una base ortogonal. Un buen ejemplo que permite acercar a las personas a la transformada de Fourier, es el oído. El oído es el órgano sensorial que permite percibir sonidos. Al oído le llegan unas ondas sonoras, la membrana del tímpano las capta y el oído interno las transforma en una descomposición de distintas frecuencias. El ser humano es capaz de discernir que frecuencias percibe a lo largo de la duración de la señal, mientras que la transformada de Fourier obtiene el conjunto de frecuencias de toda la señal. Por tanto, la transformada de Fourier no es capaz de concretar en que momento se produce una frecuencia en concreto. -3 -2 -1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 (a) Función f(x) -3 -2 -1 0 1 2 3 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 f(x) F1(X) F2(X) F3(X) F4(X) F5(X) F6(X) (b) Función f(x) -3 -2 -1 0 1 2 3 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 f(x) F1(X) F2(X) F3(X) F4(X) F5(X) F6(X) (c) Función f(x) Figura 2.1 (a)Se tiene una señal f(x) . (b)Se descompone la señal en distintas frecuencias de senos y cosenos. Esto es lo que hace el oído interno. (c)Se realiza una aproximación de la función f(x) teniendo en cuenta cuando se aplican las frecuencias y magnitud tienen. La transformada de Fourier, solo es capaz de obtener la descomposición frecuencial de la función f(x). 2.2.1 Series y Transformadas de Fourier Antes de comenzar viendo la aplicación directa que tiene Fourier sobre este trabajo, la transformación de los datos a unos ejes universales, es necesario introducir algunas propiedades y peculiaridades que tienen las series y transformadas de Fourier. Se viene hablando que las series y transformadas de Fourier tienen un carácter geométrico bastante importante que justifica su amplio uso. Esto es debido a que están relacionadas con el espacio de Hilbert. Decimos que Hes un espacio de Hilbert cuando cualquier sucesión de Cauchy converge a un elemento en el espacio, o lo que es lo mismo, que la norma de las diferencias tiende a cero. Una sucesión de Cauchy, es aquella que para una distancia determinada llega a converger en un término tal que la distancia entre dos términos cualesquiera posteriores es menos a la dada. En otras palabras, un espacio Hilbert es un espacio vectorial en el que se pueden medir la distancia entre vectores y que además es completo, lo que implica que toda sucesión de Cauchy(sucesión cuyos términos se van aproximan cada vez más entre sí) converge a un elemento dentro del propio espacio. La estructura adicional que distingue a un espacio de Hilbert es su producto interno, operación que toma dos vectores y le asigna un número. A partir de este producto interno, se puede definir una norma los vectores del espacio, pudiéndose medir distancia entre ellos. Por ejemplo, al tener el espacio euclídeo producto interno que puede expresarse como norma uno, es espacio Hilbert. El hecho de que todos los espacios finito-dimensionales con producto interno sean espacios Hilbert, permite extrapolar estos resultados a espacios infinito-dimensionales, como son las series y transformadas de Fourier o las teorías de wavelets. Para poder ver el carácter geométrico de las series y transformadas de Fourier y la importancia de que sean en espacios de Hilbert, es necesario primero ver el producto interno entre dos funciones. Sean dos funciones f(x)yg(x), se define el producto interno en el intervalo x∈[a,b]como: ⟨f(x),g(x)⟩=Zb a f(x)g(x)∗dx (2.2) donde g(x)∗es el conjugado de g(x). 2.2 Fourier y Wavelet Transforms 5 Si en vez de tomar la ecuación (2.2) como continua y se aproxima como un sumatorio acotado con los primeros kelementos, resultaría que se llega a la aproximación de Riemann siguiente: b−a n−1⟨f,g⟩= n ∑ k=1 f(xk)g(xk)∗∆x(2.3) donde n es número de puntos tomados en la aproximación. Notar que si se hace que n→∞ , o lo que es lo mismo, ∆x→0, se llega fácilmente a la ecuación (2.2). Esto conduce a poder relacionar el producto interior con el espacio de Hilbert: ∥x∥2=p⟨f,f⟩= (Zb a f(x)f(x)∗dx)1 2.(2.4) Es este producto escalar el que nos permite proyectar una función en un nuevo sistema de coordenadas de base ortogonal. Al ser espacio de Hilbert se puede llevar este resultado a espacio de infinitas dimensiones, es decir, aplicar esto a las series y transformadas de Fourier. En concreto, las series de Fourier son una proyección a una base ortogonal de las funciones seno y coseno. En ese supuesto, la función f(x)puede ser reescrita de la siguiente manera: f(x) = a0 2+ ∞ ∑ k=1 (akcos(kx)+bksin(kx)).(2.5) Para calcular los coeficientes akybkse supondrá un periodo de la función f(x)de 2π: ak=1 πZπ −π f(x)cos(kx)dx.(2.6) bk=1 πZπ −π f(x)sin(kx)dx.(2.7) En las ecuaciones (2.6) y (2.7) se puede vislumbrar que lo que se está haciendo es el producto interno entre la función f(x) y las funciones cos(kx) y sen(kx) , respectivamente. Para una mejor claridad, las ecuaciones anteriores se pueden expresar de la siguiente manera: ak=⟨f(x),cos(kx)⟩ ∥cos(kx)∥2.(2.8) bk=⟨f(x),sin(kx)⟩ ∥sin(kx)∥2.(2.9) En caso de que se quisiera expresar las series de Fourier en su forma compleja, solo habría que recordar la expresión de la fórmula de Euler: eikx =cos(kx)+ isin(kx).(2.10) Introduciendo la ecuación de Euler, las series de Fourier expresada en su forma compleja, quedaría: f(x) = ∞ ∑ k=−∞ ckeikx = ∞ ∑ k=−∞ (αk+iβk)(cos(kx) +isin(kx)).(2.11) La transformada de Fourier resulta una extensión de las series de Fourier, en el sentido de que la transformada de Fourier es aplicable a funciones genéricas no periódicas. Esto resulta de gran importancia porque en la realidad la mayoría de señales que se tratan, no tienen un carácter completamente periódico. Se verá a continuación como se llega de las series de Fourier a la transformada de Fourier. Se define el dominio de las series de Fourier, x∈[−L,L) , teniendo en mente que L→∞ , se obtendría: f(x) = a0 2+ ∞ ∑ k=1 (akcos(kxπ L)+bksin(kxπ L)).(2.12) 6Capítulo 2. Compressed Sensing Si se expresa la ecuación anterior en función de los coeficientes ck: ck=1 2L+ZL −L f(x)e−ikπx/Ldx →f(x) = ∞ ∑ k=−∞ ckeikπx/L.(2.13) Al representar f(x) de esta manera, se puede ver que el conjunto de frecuencias es kπ/L . Este conjunto de frecuencias es continuo, en tanto que L→∞. Si se denomina ω=kπ/Ly∆ω=π/L, se llega a: f(x) = l´ ım x→∆ω ∞ ∑ k=−∞ ∆ω 2πZπ/∆ω −π/∆ω f(ξ)e−ik∆ωξ dξeik∆ωx.(2.14) Tras tomar el límite de la ecuación (2.14), la transformada de Fourier de la función f(x)es Rπ/∆ω −π/∆ωf(ξ)e−ik∆ωξ dξ , que se denominará F(f(x)) . En realidad, se puede expresar como una integral de Riemann debido al término ∆ω, dando lugar a lo siguiente: f(x) = F−1(f(ω)) = 1 2πZ∞ −∞ f(ω)eiωxdx.(2.15) f(ω) = F(f(x)) = Z∞ −∞ f(x)e−iωxdx.(2.16) Estas dos integrales son conocidas como par de transformada de Fourier. La transformada de Fourier es sumamente útil para multitud de aplicaciones debido a las propiedades que esta tiene. Entre ellas se puede destacar las propiedades de linealidad, derivación en tiempo y en frecuencia, la convolución o la traslación en tiempo y en frecuencia [ 13 ]. Todas estas propiedades sirven de base para la resolución para problemas de ecuaciones de derivadas parciales, ya que los transforman en problemas de sistemas de ecuaciones ordinarios [14]. 2.2.2 Fast Fourier Transform (FFT) Hasta ahora se ha empleado las series y transformada de Fourier para funciones continuas. En la realidad, no se tendrán estas funciones bien definidas, sino paquetes discretos de información encapsulados en vectores de datos. Por tanto, se tendrá que adaptar el concepto de transformada de Fourier a transformada discreta de Fourier (DFT). La idea fundamental de este método es discretizar la función f(x) en espacios regulares, de tal manera que se le aplicará la transformada a cada una de esas discretizaciones. Sin embargo, esto en la práctica no se utiliza. Supongamos que tenemos una matriz de dimensiones n×n , el cálculo de operaciones tendrá un coste de nn . Como esto computacionalmente no es rentable, lo que se usa es la transformada rápida de Fourier (FFT), cuyo coste de computación de operaciones es nlog(n) . La idea que permite esta gran optimización, es descomponer la transformada que se quiere calcular, en otras más simples y así sucesivamente hasta legar a transformadas de dos elementos. Una vez resultas estas transformadas más simples, se han de agrupar en otras de nivel superior que han de resolverse nuevamente y así sucesivamente hasta llegar al nivel más alto. Una vez hecho esto, se procede a reordenar los resultados obtenidos. La transformada rápida de Fourier es ampliamente usada en la práctica, en la compresión en tiempo real de imagen en formato JPG o audio en formato MP3 [15], por mostrar algunos ejemplos. Para llegar al resultado matemático de la transformada rápida de Fourier (FFT), se ha de partir primero de la definición de la transformada discreta de Fourier (DFT): F(fk) = n−1 ∑ j=0 fje−i2πjk/n,(2.17) donde fjrepresenta la función fkdiscretizada en n pasos. Para un número determinado de discretizaciones n , se tiene que la frecuencia fundamental es ωn= 2.2 Fourier y Wavelet Transforms 7 e−2πi/n. Expresando la ecuación anterior en forma matricial quedaría:        F(f1) F(f2) F(f3) . . . F(fn)        | {z } F(f) =         1 1 1 ··· 1 1ωnω2 n··· ωn−1 n 1ω2 nω4 n··· ω2n−1 n . . .. . .. . ..... . . 1ωn−1 nω2n−1 n··· ω(n−1)2 n         | {z } F        f1 f2 f3 . . . fn        |{z} f .(2.18) El vector columna F(f) contiene todos los coeficientes de Fourier de la función discretizada f . Para pasar de la función entrada f a F(f) , se emplea la matriz F , que es una matriz unitaria de Vandermonde[ 16 ], que contiene elementos reales e imaginarios, es decir, magnitud y fase. Por tanto, la ecuación (2.18) puede ser reescrita como F(f) = F f . A continuación, lo que se verá ahora es la relación directa entre la DFT y la FFT. Para un mayor grado de entendimiento, se supone que se van a tratar 2048 datos, por lo que n=2048 =211 . Para empezar, si se implementara la DFT se tendría un coste de operación de n2=222 , algo inconmensurable, mientras que si se utiliza la FFT, el coste de operación es nlog(n) = 211 log(211) =≈213 . En este caso, la matriz F2048 , puede ser reescrita de la siguiente manera: F(f) = F2048 f=I1024 −D1024 I1024 −D1024F1024 0 0F1024fpar fimpar,(2.19) donde fpar y fimpar son los índices de los elementos pares e impares de f , respectivamente. Por otro lado, la matriz I1024 es la matriz identidad con dimensión 1024 ×1024 , y la matriz D1024 se expresa como: D1024 =        1 0 0 ··· 0 0ω0··· 0 0 0 ω2··· 0 . . .. . .. . ..... . . 0 0 0 ··· ω1024        .(2.20) Se puede repetir el proceso de descomponer las matrices n=2p . En caso de que los datos no se puedan expresar en potencias de dos, lo que se hace es rellenar el vector de datos con ceros hasta convertirlo de tal manera. A continuación, se mostrará una aplicación directa de la transformada rápida de Fourier (FFT), filtrado de ruido. Como el propio nombre indica, lo que se hará, será tomar una señal ruidosa, compuesta por la señal que interesa más el ruido, y transformarla de tal manera que tras distintas operaciones, aparezca la señal interesada. El motivo por el que se usa Fourier, es porque tras aplicar la FFT, se puede observar que hay elementos que pueden ser suprimidos, quedando solo con los elementos importantes. Se tiene la siguiente función continua: f(t) = sin(2πf1t)+sin(2πf2t),(2.21) donde f1=60Hz y f2=100Hz , son las frecuencias principales de la señal. A esta función continua le añadimos ruido blanco, con el fin de que parezca que es una señal ruidosa. El ruido blanco no es más que un ruido aleatorio que se caracteriza por el hecho de que sus valores de señal, no guardan una correlación estadística. El ruido blanco, además, contiene todas las frecuencias y muestran la misma potencia. Se pretende, ahora, ver la efectividad de filtrado. Se sabe que tras transformar, se aplique el filtro y se haga la inversa de la transformada, el resultado debe ser algo similar a la función f(t) . Conviene aclarar que se hace en la parte del filtro. Lo que se pretende es ver la densidad de potencia espectral (PSD), que viene de normalizar las componentes de Fourier y ver que potencia corresponde a qué frecuencia, con el fin de escoger un cierto valor a partir del cual, poder cortar aquellas frecuencias que no interesan. Los resultados son los siguientes: 8Capítulo 2. Compressed Sensing 2.2.3 Wavelets Las wavelets [ 17 ] son señales, o formas de onda, que tienen una duración temporal limitada y media cero. Para poder ver la motivación detrás de este concepto, es necesario primero mencionar el principio de incertidumbre. Es importante ya que limita la capacidad que se tiene para medir a gran resolución tiempo y frecuencia simultáneamente. Por ejemplo, una función continua en el tiempo, permite obtener información del tiempo en todo su dominio, sin informar nada de la frecuencia. El caso inverso sucede con Fourier, se tiene todo el contenido frecuencial, pero no el tiempo en el que ocurren esas frecuencias. Esta incertidumbre lleva a buscar una solución intermedia, en la que se obtenga un contenido de calidad tanto en frecuencia como en tiempo. Eso es lo que otorgan las wavelets. Las wavelets se introducen debido a que permiten hacer el análisis de Fourier a bases ortogonales más generales y, más importante, son capaces de manejar de manera destacable el principio de incertidumbre. Las wavelets se utilizan en distintos campos como la neurociencia [ 18 ] o en los fenómenos turbulentos [ 19 ]. La idea básica detrás de todo esto, es comenzar con una función conocida ψ(t) , llamada función wavelet 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 Tiempo [s] -15 -10 -5 0 5 10 Amplitud Con ruido Sin ruido 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 Frecuencia [Hz] -2 -1 0 1 2 PSD Sin ruido Filtrada 0 50 100 150 200 250 300 350 400 450 500 Tiempo [s] 0 100 200 300 Amplitud Threshold Con ruido Filtrada Figura 2.2 Filtrado de ruido con FFT. (arriba)La señal limpia sin ruido y la señal ruidosa, compuesta por la señal sin ruido más cierto ruido aleatorio.(medio) El filtrado se hace de tal manera que nos quedamos con las frecuencias que tengan una energía mayor que un umbral definido previamente.(abajo) Una vez filtrado y hecha la transformada inversa, comparando resultados. 2.2 Fourier y Wavelet Transforms 9 madre, a partir de la cual se genera una familia de funciones de condiciones parecidas a esta. Esta función es de la forma: ψa,b(t) = 1 √aψ(t−b a),(2.22) Los parámetros a y b son los encargados de relacionar la wavelet madre con la familia de funciones. Para mostrar en que consiste este proceso, se mostrará el ejemplo más sencillo de wavelet. En 1909 Alfred Haar, introdujo lo que hoy en día se conoce como la primera wavelet, a pesar de que en tal fecha no existía aún ese concepto. La función wavelet madre de las funciones de Haar ψ(t) puede ser descrita como: ψ(t) =    1si 0≤t<1/2, −1si 1/2≤t<1, 0si no (2.23) A continuación, se presenta un ejemplo gráfico de tres wavelet de Haar: -0.2 0 0.2 0.4 0.6 0.8 1 1.2 t -1.5 -1 -0.5 0 0.5 1 1.5 1,0 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 t -1.5 -1 -0.5 0 0.5 1 1.5 0.5,0 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 t -1.5 -1 -0.5 0 0.5 1 1.5 0.5,0.5 Figura 2.3 Tres wavelets de Haar para los primeros dos niveles. 10 Capítulo 2. Compressed Sensing 2.3 Compression Lo visto hasta ahora, Fourier y wavelets, son métodos que sirven para generalizar dimensiones de espacio elevadas a dimensiones de espacio pequeñas, como puede ser 2D o 3D. Tanto Fourier como las wavelets han servido para proveer de algoritmos robustos en el proceso de filtrado y compresión de imagen, lo que permite una buena maniobrabilidad de grandes dimensiones de datos, algo que, hasta el momento, era inviable computacionalmente. A pesar de que la mayoría de compresión ha sido utilizada para imágenes y audio, también tiene gran aplicación en el campo de la ingeniería. Algunos de estos ejemplos pueden ser la resolución de los sistemas de ecuaciones diferenciales parciales, que pueden aparecer en la resolución de problemas físicos como distribución de calor [ 20 ], y también se usa para caracterizar de manera instantánea, con una representación sparse, los fenómenos turbulentos, que es donde se adentrará este trabajo. Antes de comenzar con como ver como se puede hacer esto último, conviene mostrar un ejemplo de como se puede comprimir una imagen, quedando solo los términos que tienen una mayor representación. Esto sirve para introducir brevemente el concepto de sparse, donde solo unos pocos son no nulos y, por tanto, son los relevantes. Anteriormente, se ha supuesto que las imágenes son compresibles sin explicar de donde viene. Si suponemos que se trata de una imagen en blanco y negro de dimensión 30 ×30 , donde los píxeles pueden tomar únicamente valores binarios, las posibilidades que se tienen de estructurar la foto de manera distinta es de 2900 , que por poner un ejemplo de tal dimensión, es mayor que el número total de partículas conocidas en nuestro universo. Además, este número se puede hacer aún mayor si se tiene en cuenta las imágenes a color y que las dimensiones con las que se tratan hoy día, son del orden de los megapíxeles. Todo lo anterior conduce a que las imágenes que puede tener en mente, representan una parte ínfima del total de las imágenes. Debido a las infinitas posibilidades de generar imágenes distintas, la mayoría de imágenes pueden considerarse aleatorias o, al menos, pertenecen al conjunto de imágenes que nosotros no percibimos, es decir, no representan nada que para nosotros tenga algún sentido. Por tanto, las altas dimensiones quedan para representar estas imágenes aleatorias. Esto último implica, que si se quiere representar una imagen que nosotros entendemos, la mayoría de dimensiones son redundantes y no necesarias. Es aquí, dónde nace la idea y se puede explotar el concepto de compresión de imagen. En el siguiente ejemplo, se va a llevar a cabo la compresión de una imagen de dos maneras distintas, mediante FFT y wavelets. Para llevar a cabo dicha compresión con Fourier, se aplica la transformada rápida de Fourier en dos dimensiones distintas. Si se representa la imagen como una matriz, la FFT se aplica en filas primero y en columnas posteriormente. El resultado que se obtendría si se hiciera primero en columnas y posteriormente en filas sería el mismo. Se mostrarán a continuación los resultados en las ilustraciones de la figura 2.4. Ahora se verá lo que resultaría al aplicar wavelets como método de compresión. Mencionar únicamente que lo que se aplica es la transformada de wavelet discreta (DWT), que se deduce igual que la transformada discreta de Fourier. Los resultados se reflejan en las ilustraciones de la figura 2.5 Se puede observar que ambas transformaciones logran una compresión eficaz, económica y precisa, permitiendo eliminar información no relevante y, por tanto, reduciendo la dimensionalidad de las imágenes. En la próxima sección se ampliará el concepto de compresión a una idea mucho más general, que consiste expresar una cantidad masiva de información en función de unos pocos de datos, llamados sparse. 2.3 Compression 11 Imagen Original (a) Compresión al 95% (b) Compresión al 99% (c) Compresión al 99.8% (d) Figura 2.4 Compresión de imagen a distintos niveles utilizando la transformada rápida de Fourier (FFT) bidimensional. Se puede ver que dejando solo el 5% de la imagen original, la imagen apenas pierde información. A medida que la compresión se hace más agresiva, los detalles de la imagen original desaparecen. 18 Capítulo 2. Compressed Sensing 2.5.2 Compresión de una imagen Se va a tener una imagen de alta calidad, la cual se pretende comprimir con el fin de que ocupe menos en memoria. Para ello lo que se hace es hacer la transformada rápida de Fourier para las filas y columnas de la imagen, se hará una umbralización con el fin de obtener los términos más representativos y realizar la transformada inversa, resultando: Imagen Original Componentes Fourier Completas Componentes Fourier Parciales Reconstrucción Figura 2.9 Se tiene el proceso de compresión de una imagen en escala de grises. (arriba derecha)Se puede ver como están dispuestos las componentes de la transformada rápida de Fourier (FFT) en escala logarítmica. Para interpretar esto, basta con entender que las imágenes representan baja frecuencia y los píxeles más claros son los de mayor importancia. (abajo izquierda)Esto implica que al tener en cuenta el umbral de componentes deseadas, la mayoría de términos se condensen en el centro de la imagen. 2.5 Ejemplos 19 2.5.3 Optimización norma uno vs. mínimos cuadrados Este último ejemplo resulta ilustrativo a la hora de entender por qué es interesante usar la norma uno frente a mínimos cuadrados en cuanto a sparse se refiere. Lo que se va a hacer es resolver el siguiente sistema y=Θs , donde Θ∈p×n con p=100 filas y n=500 columnas. En general existen infinitas soluciones que resuelven este sistema. Sin embargo, existe una pequeña posibilidad en la que las filas de la matriz Θ sean linealmente dependientes mientras las columnas no sean consistentes con esas filas. Esta pequeña posibilidad puede disminuirse si se aumentara el número de columnas, ya que generaría más soluciones. Para resolver este problema indeterminado se usará minimización norma uno y minimización norma dos: 0 100 200 300 400 500 -0.2 -0.1 0 0.1 0.2 si Norma L1 0 100 200 300 400 500 -0.2 -0.1 0 0.1 0.2 Norma L2 -0.1 -0.05 0 0.05 0.1 0 100 200 300 400 500 600 Histograma(si) -0.1 -0.05 0 0.05 0.1 0 50 100 150 200 250 300 350 400 Figura 2.10 Comparación entre resolución mediante minimización norma uno y minimización norma dos. (izquierda) Se puede observar un único valor que es predominante y el resto es de valor muy bajo, casi cero. Esto lo que hace es impulsar que el vector ssea sparse. (derecha)La minimización norma dos se puede ver que no es tan buena, ya que la mayoría de términos son de valores parecidos, sin predominar unos pocos. Es por esto mismo por lo que conviene usar norma uno.. 3 Optimal Sparse Sensing La estrategia expuesta del compressed sensing, es conveniente cuando la reconstrucción de una señal tiene la dimensión de datos es elevada y cuya estructura es desconocida de antemano. Por ejemplo, si se quisiera hacer la reconstrucción de una imagen de una playa, de una cara y de una casa, compressed sensing resulta ideal, ya que ofrece una solución bastante aceptable. Por tanto, esta estrategia resulta flexible para cuando tenemos una señal de un contenido incierto. En cambio, si se sabe que la señal va a ser únicamente imágenes de una casa, es posible hacer el diseño de unos sensores óptimos, teniendo una base de transformación ajustada a ese tipo de datos en concreto. Si se tienen unos datos ajustados, se pueden extraer las características dominantes de ellos, como lo hace la descomposición ortogonal adecuada, en inglés Proper Orthogonal Descomposition (POD) [ 24 ]. Estas características de rango bajo, que pueden ser inferidas de los datos de la señal, facilitan el diseño óptimo de los sensores, siendo estos los adecuados para el problema que se tiene entre manos. El éxito de la descomposición ortogonal adecuada (POD) se debe a la existencia de comportamientos significativos que están en patrones de dimensión reducida. Esta técnica lo que busca es aprovechar esto, para producir sistemas dinámicos de orden inferior que sean capaces de modelar de manera precisa el sistema dinámico al completo. La POD sirve de apoyo al campo de modelos de orden reducidos (ROM), ya que proyecta el problema de ecuaciones diferenciales parciales en subespacios de orden reducido, donde estas ecuaciones pueden ser resueltas con mayor facilidad. A su vez, los llamados métodos de interpolación empírica (EIM), requieren que los datos estén representados en una base que tenga características POD. Lo que se verá en este capítulo será toda la formulación matemática que se requiere para abordar el problema de la optimización de sensores, como la brevemente comentada descomposición ortogonal adecuada (POD), pero también la búsqueda de los distintos algoritmos que nos permiten encontrar la solución del problema. Por último, se concluirá con unos ejemplos de aplicación de diversas índoles, todos con el mismo concepto en mente, como es el uso de sensores óptimos para la reconstrucción. 3.1 Proper Orthogonal Descomposition (POD) La descomposición ortogonal adecuada (POD), es una de las técnicas de reducción de dimensionalidad más importantes que se conocen hoy en día. POD lo que hace es expresar un dato de grandes dimensiones x∈Rn , como puede ser una imagen o el flujo de un fluido, como una combinación lineal de ciertos modos propios ortogonales, que definen un espacio de dimensión baja. Los estados se proyectan en estos POD subespacios, lo que permite un mejor tratamiento de datos, algo que sería impensable si se quisiera trabajar con el espacio completo. Sin embargo, esto no se puede hacer de manera instantánea, ya que es necesario tener una base de datos de entrenamiento extensa con la que ajustar la base POD al problema que se requiera abordar. A modo de ejemplo, si se quisiera reconstruir una cara de una persona con pocos datos, es necesario tener muchas imágenes de una persona como datos de entrenamiento. Lo que se obtendría en este ejemplo, sería un subespacio dónde los modos del POD son la nariz, los ojos o la boca. Para realizar una representación de dimensión baja de x en una base POD, lo que se busca es expresar la siguiente relación lineal: xi≈ r ∑ k=1 ak(ti)ψk(x),(3.1) 21 22 Capítulo 3. Optimal Sparse Sensing donde los coeficientes ak(ti) dependen del tiempo y los términos ψk(x) solo dependen del espacio, por lo que representa una descomposición espacio-tiempo. Estos términos pueden obtenerse fácilmente al aplicar la descomposición singular de valores (SVD), que se explicará detenidamente en la subsección 3.1.1. A modo resumen, lo que se pretende es expresar Xcomo: X=ΨΣVT≈ΨrΣrVT r,(3.2) donde las matrices Ψr y Vr contienen las primeras rcolumnas del vector singular izquierdo y derecho respectivamente, y la matriz Σres la matriz identidad cuyos elementos son los valores singulares. Como se busca reducir la dimensión del estado X , lo que se intenta es expresar la siguiente proyección ortogonal a=ΨT rx . Por tanto, se puede ver que expresándolo en una base POD, el problema se simplifica en cuanto a dimensión de manera considerable, lo que permite llegar a la solución en multitud de problemas. Lo que se tendrá que ver ahora es cómo se puede llegar a esa descomposición de valores singulares y la elección óptima del valor r. Aunque en la práctica lo que se hace es tomar aquellos primeros rvalores que captan un cierto porcentaje de los datos totales, también se verá la elección de un umbral más preciso. 3.1.1 Singular Value Descomposition (SVD) La descomposición en valores singulares, en inglés Singular Value Descomposition (SVD), es una de las factorizaciones matriciales más importante de los últimos tiempos. Los usos principales de esta factorización son expresar una matriz en otras que tengan un rango inferior, con el fin de simplificar el problema, y calcular la pseudo-inversa en el caso de que se tengan matrices no cuadradas, siendo esto importante cuando el sistema Ax =bno es determinado. Anteriormente, se ha explicado de manera detenida la transformada rápida de Fourier, en la sección 2.2.2. El SVD no es más que una generalización de la FFT. Generalización en el sentido de que la transformada rápida de Fourier se acomete en situaciones idealizadas, mientras que la descomposición en valores singulares es aplicable en la mayoría de casos de manera correcta. A su vez, en lo que al tratamiento de datos se refiere, la FFT proporciona una base genérica, mientras que el SVD da una base ajustada, una base POD. Es por esto último, por lo que se usará en este capítulo. En los problemas en los que los datos son multidimensionales, el SVD permite encontrar los patrones relevantes de una dimensión mucha menor. Esto es algo muy potente, ya que no es necesario que un experto de la materia en concreto tenga que observar los datos e inducir algo a partir de ellos, sino que simplemente, con un tratamiento puramente numérico de los datos, es posible hacer lo mismo. Además, este tratamiento es numéricamente estable y aplicable a cualquier matriz. La reducción del orden es una aplicación fundamental del SVD, sin embargo, no es la única. El SVD también es usado para obtener la pseudo-inversa de matrices no cuadradas. Esto permite resolver sistemas Ax =b que sean incompatibles y compatibles indeterminados. A su vez, esta resolución puede entenderse como una transformación geométrica de autovectores y autovalores, en los que los ángulos y distancias se preservan. Se tiene una matriz de datos X∈Cn×m:     . . .. . .. . . x1x2··· xm . . .. . .. . .    .(3.3) Cada columna xk∈Cn , representa un experimento o simulación. Por ejemplo, en el caso de que se esté midiendo la vorticidad de un fluido, cada columna recoge fotografías de la vorticidad en distintos periodos de tiempo. Algo similar se hace con las imágenes. Si se tienen imágenes de distintas personas de dimensión 1024×1024 , cada columna es una imagen de una persona de tamaño 1048576×1 . En este sentido, se tienen mdistintos experimentos. Cada columna es normalmente conocida como instantánea, donde mes el número de instantáneas que se tienen del experimento en X. 3.1 Proper Orthogonal Descomposition (POD) 23 El SVD transforma la matriz X∈Cn×men las siguientes matrices: X=UΣV−1,(3.4) donde U∈Cn×n y V∈Cm×m son matrices unitarias y ortonormales, mientras que Σ∈Rn×m es una matriz real de tal forma que en la diagonal se encuentran valores no negativos y nulos debajo de la misma. En muchos casos se tendrá que la matriz de datos X es alta y delgada, o lo que es lo mismo, se tienen datos extensos pero pocas instantáneas de los mismos, n≥m . Esto implica que la matriz Σ tienen mvalores no nulos en la diagonal, por lo que esta matriz puede ser reescrita como Σ=ˆ Σ 0. Si esto último se introduce en la ecuación (3.4), se obtiene lo que se conoce como SVD economizado: X=UΣV−1=ˆ Uˆ UTˆ Σ 0V−1=ˆ Uˆ ΣV−1(3.5) El nombre de descomposición en valores singulares es debido a que la matriz ˆ Σ contiene en su diagonal los valores singulares. A su vez, las columnas U y V son conocidas como los vectores singulares izquierdo y derecho, respectivamente. Es conveniente usar el SVD economizado, debido principalmente a que las operaciones que se tienen que realizar son menos costosas, computacionalmente. De no usarse, se estaría utilizando matrices de orden muy alto, mientras que en su interior la información relevante tiene una dimensión mucho menor, por lo que se está despreciando coste computacional. Sin embargo, el uso del SVD no llega a ser suficiente. Esto es debido a que existen valores singulares en ˆ Σ que tienen un mayor peso que otros. Esto motiva que una vez obtenido el SVD economizado se trunquen ciertos valores singulares según la energía que tengan sobre el total. Para truncar hay que tener en cuenta muchos factores, entre los que hay que destacar las especificaciones que se quieran cumplir, el ruido que pueda existir y la distribución de los valores singulares. Existen muchas formas de realizar esto, tanto truncamientos simples como complejos. Lo que se suele hacer es truncar de tal manera que se capture un nivel energía de los datos originales, un cierto porcentaje. Aunque parezca poco laborioso y simple, es muy usado este método. Sin embargo, se han realizado distintos trabajos para obtener un truncamiento óptimo. Entre todos estos, destacar el trabajo realizado por Gavish y Donoho [ 25 ]. En dicho trabajo, lo que se calcula es el umbral óptimo a partir del cual truncar. Para ilustrar todo lo comentado anteriormente, se propone el siguiente ejemplo. Lo que se hace en este en las ilustraciones de la figura 3.1 es descomponer una imagen haciendo uso del SVD optimizado, ver una representación de los valores singulares y las dos posibilidades de truncamiento. El SVD sirve de base para muchos procedimientos donde se realiza una reducción de dimensionalidad. Estos métodos abarcan desde el análisis de componentes principales, en inglés Robust Component Analysis (PCA), usado en el campo de la estadística [ 26 ], funciones ortogonales empíricas usadas en climatología [ 27 ], la descomposición ortogonal adecuada en dinámica de fluidos [ 28 ] y el análisis de corrección canónica, en inglés Canonical Correlation Analysis (CCA), [29]. Es posible que el uso más potente que se le puede dar a una de las propiedades del SVD es que proporciona una aproximación óptima de rango inferior a una matriz X . De hecho, el SVD proporciona aproximaciones de bajo rango de manera jerárquica, de tal modo que una aproximación de rango rse obtiene al mantener los principales rvalores y vectores singulares, descartando el resto. Gracias a la ecuación (3.5) , donde se establece una notación de una base truncada de SVD, se va a usar el hecho de que Σ es una matriz diagonal, por lo que la aproximación del SVD de rango rse obtiene mediante la suma de rmatrices distintas de rango unidad: ˆ X=ˆ Uˆ ΣV−1= r ∑ k=1 σkukv−1 k=σ1u1v−1 1+σ2u2v−1 2+···+σrurv−1 r.(3.6) Para un rango rdado, no hay una mejor aproximación para X , en el sentido de norma dos (mínimos cuadrados), que la aproximación truncada de SVD ˆ X . Por lo tanto, los datos de alta dimensionalidad pueden describirse mediante unos pocos patrones dominantes dados por las columnas de ˆ Uyˆ V. Esta propiedad del SVD es fundamental. Existen numerosos ejemplos de conjuntos de datos que contienen mediciones de alta dimensionalidad, lo que resulta en una matriz de datos X de gran tamaño. Sin embargo, a menudo hay patrones dominantes de baja dimensionalidad en los datos, y la base truncada de SVD ˆ U , aunque también se usará Ψr , proporciona una transformación de coordenadas desde el espacio 24 Capítulo 3. Optimal Sparse Sensing 0 1000 2000 3000 4000 r 100 101 102 103 104 105 106 Valor singular, r Valores singulares con truncamiento al 90% 0 1000 2000 3000 4000 r 100 101 102 103 104 105 106 Valores singulares con truncamiento óptimo Figura 3.1 Representación valores singulares de una imagen. A izquierda, se puede observar el truncamiento cogiendo el 90% de energía de los valores singulares. A derecha, la otra opción es truncamiento óptimo, usando umbral =4/√3√nσ con nnúmero de filas y σ cero, indicando que se desconoce el ruido. Por debajo de ese umbral, se desecharán los valores singulares. de medición de alta dimensionalidad hacia un espacio de patrones de baja dimensionalidad. Esto tiene la ventaja de reducir el tamaño y la dimensión de los conjuntos de datos grandes, lo que proporciona una base manejable para el análisis. 3.2 Sensor Placement for Reconstruction Se ha visto que es posible reducir la dimensionalidad del problema al tratar los datos, tras disponerlos en una base determinada, en concreto, en una base POD. La idea que puede aparecer tras saber esto, es si es posible obtener de manera directa los datos reducidos en dimensión, con lo que se ahorraría un paso importante computacional. Siguiendo esta filosofía, si ponemos el ejemplo la vorticidad de un fluido a lo largo del tiempo, interesa saber si es posible determinar la ubicación óptima de sensores con los que reconstruir por completo la vorticidad del fluido. De primeras, se tendría que comprobar caso por caso si la ubicación del sensor es la adecuada, lo que conlleva un problema combinatorio de encontrar xsensores en las nposibilidades. Esto resultaría algo intratable computacionalmente. Sin embargo, como se verá en esta sección, será posible abordar el problema. Tanto las situaciones en las que los sensores son pocos y costosos, como en las que los sensores son muchos y baratos, aprovecharán esta solución del problema. Esta será la aplicación directa que se le dará a toda esta matemática explicada previamente, sin embargo, puede tener multitud de usos, como por ejemplo utilizar los conceptos sparse para llevar a cabo una clasificación [30]. Se persigue optimizar la ubicación de los sensores para reconstruir todo el estado de alta dimensión gracias a las medidas tomadas en esos sitios concretos. Para la reconstrucción es necesario conocer en qué base ajustada se expresan los datos, ya que la reconstrucción del estado de alta dimensión se expresa: x=Ψra.(3.7) Ψr∈Rn×rrepresenta la base ajustada, como puede ser una base POD que se puede obtener fácilmente tras aplicar el SVD, a∈Rr son los coeficientes de la base ajustada y x∈Rn es el estado de alta dimensión. Si la 3.2 Sensor Placement for Reconstruction 25 anterior ecuación (3.7) es expresada en forma coordenada, se obtiene: xj= r ∑ k=1 Ψjkak,(3.8) donde Ψjk es una matriz que expresa en sistema de coordenadas la base ajustada Ψr . Lo que se busca es encontrar de manera efectiva una matriz de medidas C∈Rp×n que está optimizada de tal manera que es capaz de recuperar los coeficientes a de las salidas de los sensores y . Estas salidas no son más que los p elementos escogidos de x y=Cx = (CΨr)a.(3.9) Expresando esta ecuación (3.9) en sistema de coordenadas: yi= n ∑ j=1 Ci jxj= n ∑ j=1 Ci j r ∑ k=1 Ψjkak.(3.10) Cuando se desconoce el estado x , es posible reconstruirlo mediante una aproximación de los coeficientes de base ajustados, a: ˆ x=Ψrˆa,donde ˆa=Θ−1y= (CΨr)−1ysi p=r Θ∗y= (CΨr)∗ysi p>r(3.11) En la ecuación (3.11) hay que destacar que Θ∗ es la pseudo-inversa, ya que la matriz no es cuadrada. Esto es debido principalmente a que existe un mayor número de sensores,p, qué modos usados en la reconstrucción, r. Se persigue buscar la ubicación óptima de los sensores que permiten la mejor reconstrucción posible del estado ˆ x . Siguiendo la misma línea, la matriz de medidas C intenta escoger las mejores filas de Ψr que permiten la mejor reconstrucción, pero también, que condicione de manera óptima la inversión dé Θ . Respecto a esta manera óptima de inversión, aparece el número de condición [ 31 ]. El número de condición de la matriz Θ mide cuanto de sensible es la matriz a las multiplicaciones e inversiones cuando hay cierto ruido en la entrada. El número de condición es el cociente entre el máximo y mínimo valores singulares de Θ: κ(Θ) = σmax(Θ) σmin(Θ).(3.12) Si el número de condición es grande, entonces el error se amplifica. Es por esto mismo, por lo que se busca que el número de condición sea lo más cercano a uno. Por tanto, este valor ha de ser estudiado detenidamente, de tal manera que la matriz de medidas C ha de ser seleccionada de tal manera que haga tender este valor a uno. Se puede observar que la estructura que debe tener C ha de escogerse de manera minuciosa para que cumpla ciertos requisitos. En lo que sigue a continuación, se supondrá que M=ΘTΘ para p>r y M=Θ para p=r . Controlando el contenido de la matriz M se controla de manera indirecta el número de condición, ya sea viendo el determinante, la traza o el radio espectral. El criterio de radio espectral de M−1 hace maximizar el valor singular más pequeño de M , con el fin de buscar un valor de número de condición, cercano a uno: γ=argmin ∥M−1∥2=argmax σmin(M)con ∥λ∥=p.(3.13) La traza debe optimizarse como: γ=argmax tr(M) = argmax ∑ i λi(M)con ∥λ∥=p.(3.14) Y por último, el determinante: γ=argmax |detM|=argmax ∏ i|λi(M)|con ∥λ∥=p.(3.15) Intentar optimizar estos problemas de manera directa, ocasionaría un problema combinatorio de n sobre p, que al igual que se ha comentado anteriormente, es algo que se debe evitar a toda costa, ya que computacionalmente es intratable. A raíz de esto han surgido muchos métodos para tener una matriz óptima 26 Capítulo 3. Optimal Sparse Sensing de M , y, por tanto, indirectamente de C . Existen muchos métodos experimentales que proporcionan solución a este problema. La primera solución que apareció fue realizar un submuestreo aleatorio [ 32 ], [ 33 ] cuyos algoritmos fueron usados para la instalación de sensores marinos [ 34 ]. Más tarde, aparecieron los EIM, Empirical Interpolation Method, [ 35 ] los discretos DEIM [ 36 ] y los Q-DEIM [ 37 ] que usan el pivotaje de la matriz Q en la factorización QR. De una manera u otra, estos algoritmos proveen un muestreo óptimo usando interpolación. En lo que sigue, se explicará uno de estos algoritmos, en concreto uno que se basa en Q-DEIM. 3.2.1 Quadrature-Discrete Empirical Interpolation Method (Q-DEIM) A pesar de que lo más sencillo es tener el mismo número de sensores que de modos, p=r , existen casos en los que puede resultar conveniente introducir un mayor número de sensores. Como se verá, es posible resolver este problema de un sistema indeterminado [ 38 ]. Teniendo los datos expresados en una base POD, solo basta con aplicar factorización QR con pivoteo de columnas [ 39 ]. La factorización QR consiste en descomponer la matriz A∈Rm×ncomo: APT=QR si p=r(3.16) donde la matriz Q es ortogonal, la matriz R es triangular superior y P es la matriz de permutación. Con esta factorización se puede resolver de manera óptima la ecuación del determinante (3.15) . La manera de resolverlo es maximizando el volumen de la submatriz, ya que el volumen de la matriz coincide con el determinante. La factorización QR con pivoteo aumenta el volumen de la submatriz construida al seleccionar cada vez una nueva columna de pivoteo de máxima norma dos y posteriormente restando de cada columna su proyección ortogonal sobre la columna de pivoteo. El pivoteo lo que hace es empujar a la submatriz a tener una estructura diagonal, lo que potencia un aumento del volumen: σ2 i=|rii|2≥ k ∑ j=i|r jk|2; 1 ≥i≥k≥(3.17) |detA|=∏ i σi=∏ i|rii|.(3.18) La factorización QR con pivoteo conduce a p=r sensores (pivotes) que muestrean de la mejor manera los rmodos de Ψr , en el caso de que se dispongan el mismo número de sensores que de modos. La factorización QR con pivoteo funciona también para el caso de sobremuestreo, más sensores que modos. Por tanto, se tendrán dos casos: (Ψr)PT=QR si p=r (ΨrΨT r)PT=QR si p>r,(3.19) donde PT coincide con lo que se viene llamando matriz de medidas CT . La única diferencia entre ambos casos es la carga computacional que requieren. El caso de que p>res más costoso, ya que se necesita factorizar matrices de n×n, mientras que en el caso p=rla factorización es para una matriz de n×r. En la mayoría de situaciones n≪r . La selección de filas dada por los primeros ppivotes hacen aumentar los primeros r valores de ΘΘT , por tanto, incrementando el determinante de M=ΘTΘ . En definitiva, este método como máximo tendrá que realizar una factorización QR del orden de n3 , dando como resultado una lista ordenada de mayor a menor de npivotes, con los primeros ppivotes que serán usados para la reconstrucción en Ψr cuando p>r. A pesar de que se esté haciendo uso de estos algoritmos para implementar la ubicación óptima de sensores para la reconstrucción, existen otras aplicaciones en las que puede resultar provechoso hacer uso de estos algoritmos [40], [41]. 3.2.2 Sensores Aleatorios En general, se pueden utilizar sensores colocados aleatoriamente para estimar los coeficientes de la base ajustada a . Sin embargo, cuando p = r y el número de sensores es igual al número de modos, el número de condición suele ser muy grande. De hecho, la matriz Θ suele ser numéricamente singular y el número de condición es muy alto. Hay que recordar que un número de condición alto lleva a inestabilidad numérica, se amplifican los errores y se encuentran dificultades para invertir matrices. El sobremuestreo, que consiste en instalar más sensores que modos, p > r, hace mejorar rápidamente el número de condición, consiguiendo, generalmente, un rendimiento de reconstrucción mucho mejor. Por tanto, es algo que hay que tener en 3.3 Ejemplos 27 cuenta, si se está dispuesto a aumentar el número de sensores para cumplir unas condiciones de error de reconstrucción preestablecidas. 3.3 Ejemplos A continuación, se verán aplicaciones de la ubicación óptima de sensores para reconstrucción en tres áreas distintas. El primer ejemplo tratará acerca de la reconstrucción facial, donde el dataset son imágenes de diferentes personas. El segundo ejemplo intentará conocer la temperatura de los mares y océanos, con la distribución de unos pocos sensores. Por último, el tercer ejemplo tratará con datos de enorme dimensión, como es la mecánica de fluidos, para disponer ciertos sensores y reconstruir el campo fluido de manera bastante acertada. En los ejemplos se hará uso principalmente del algoritmo Q-DEIM, subsección 3.2.1, que en teoría es el que mejor funciona cuando se tiene una base ajustada. Sin embargo, para comparar resultados, se va a hacer uso del Compressed Sensing. En este sentido, se va a usar una metodología muy similar a la empleada en el ejemplo del capítulo anterior 2.5.1, usando optimización con norma uno como optimización con norma dos, mínimos cuadrados. 3.3.1 Extended Yale B Eingenfaces El tratamiento de imagen requiere el uso de grandes dimensiones de datos. Las imágenes tomadas por los distintos dispositivos fotográficos contienen pequeñas partes de información llamadas píxeles, que, puestas en conjunto, conforman la imagen. Hoy en día es posible encontrar dispositivos capaces de capturar una imagen de 100 megapíxeles. A pesar de esta vastedad, la mayoría de la información que contienen estos píxeles puede ser desechada. Es por esto mismo, por lo que las imágenes puedes ser comprimidas [42]. La base de datos que se empleará en este ejemplo es Yale B face [ 43 ]. Esta base de datos es de gran utilidad cuando se quiere probar reconocimiento facial y sobre todo, sirve para mostrar la tremenda reducción de dimensión que es posible llevar a cabo en imágenes. En cuanto a la estructura del set de datos, contiene 64 fotografías de una persona con distinta iluminación cada captura y un total de 38 personas distintas. Para el entrenamiento, se tomarán 32 personas de forma aleatoria, dejando el resto para el set de comprobación. Este entrenamiento será necesario para obtener la base ajustada, base POD. Para empezar, es necesario obtener la reducción de dimensionalidad. Para ello se aplica el SVD economizado y un posterior truncamiento óptimo. En la ilustración de la figura 3.2, se realiza el truncamiento óptimo, donde se hace uso de la fórmula proporcionada por Donoho [ 25 ]. Este valor singular, óptimo, r=166 indica el rango intrínseco que tiene el set de entrenamiento. A partir del valor singular óptimo, se encuentran los valores que no aportan información relevante, solo añaden complejidad dimensional al problema. 0 100 200 300 400 500 600 r 10-2 10-1 100 101 102 r Valores Singulares Modo 166 Modo 400 Modo 50 Modo 100 Figura 3.2 Representación valores singulares del set de entrenamiento. El truncamiento óptimo llega con el modo 166, a partir del cual, los valores singulares carecen de información importante. 34 Capítulo 3. Optimal Sparse Sensing Ejemplo simulación aleatoria -1012345678 2 1 0 -1 -2 Reconstrucción con 42 sensores aleatorios -1012345678 2 1 0 -1 -2 Reconstrucción con 42 sensores QR -1012345678 2 1 0 -1 -2 Figura 3.10 Reconstrucción con 42 sensores. La reconstrucción llevada ahora a cabo, resulta más idónea. A pesar de que la distribución de los sensores aleatorios no se acercan, en absoluto, a los sensores QR, ambos tienen un error de reconstrucción semejante a la figura 3.8. 4 Estimación del Viento en Zonas de Inspección Co mo se ha podido observar anteriormente, llevar a cabo la implementación de este conjunto de algoritmos puede resultar ventajoso en ciertos escenarios. Estos escenarios abarcan desde llevar la optimización de la ubicación de sensores a la agricultura [ 48 ] hasta el campo de la sismología [ 49 ]. Este trabajo pretende ir un paso más allá, respecto las aplicaciones que actualmente se han puesto en práctica. Cuando se quiere llevar a cabo una inspección de la estructura de una fábrica o de una industria, se requiere de un utillaje específico. Hoy en día, la tendencia es hacer que estas revisiones se hagan de manera remota. Algunas de las ventajas que tiene realizar este trabajo de esta manera son el ahorro de tiempo, de costos y reducción de riesgos para las personas que están haciendo la inspección. Para llevar a cabo esta labor, se está apostando por el uso de vehículos aéreos no tripulados o drones, más concretamente multirrotores. El viento puede tener un impacto significativo en la estabilidad y el control del dron. Fuertes ráfagas de viento pueden empujar al dron fuera de su trayectoria prevista y también provoca mayor resistencia de vuelo, lo que agota la batería rápidamente. La protección contra el viento y saber como actuar frente a él, es una de las tareas pendientes hoy en día. Centrando la atención en el viento, las distintas predicciones de meteorología que pueden ser conocidas, dan datos importantes, pero no decisivos, a la hora de decidir si es posible realizar un vuelo seguro con un dron. Esto es debido, a que las predicciones tienen un carácter general, no van al detalle de determinar que rango de viento va a hacer en la inspección que se quiera llevar a cabo. Por tanto, puede resultar interesante, ser capaz de reconstruir el flujo de aire del emplazamiento en cuestión, con el fin de tener una predicción precisa sobre el terreno que se va a realizar la inspección. Una vez se tiene esta estimación, es posible usarla para la planificación del camino que tendrá el UAV inspector, o bien, sirve de entrada al sistema de control con el fin de corregir las posibles desviaciones que pueda provocar el viento. A continuación, se hará una serie de recapitulaciones para la implementación de estos procedimientos. Se parte de que el terreno de inspección puede ser conocido de antemano, pudiendo llevarse a cabo la generación de distintas simulaciones, con cierto programa, para distintos casos de velocidad y orientación de viento. A estas simulaciones se les aplica los algoritmos de reducción de dimensionalidad y obtención óptima de sensores. La solución que se va a dar una vez se tenga resuelto el problema, son el número y la ubicación de sensores necesarios, para realizar la reconstrucción. La relación que tiene esto con un escenario global es la siguiente: cada ubicación de sensor que salga en la solución va a servir para colocar en ese lugar un pequeño UAV que tenga implementando un sistema de medición de velocidad y orientación del viento. Una vez estos drones estén en disposición, tomarán las mediciones pertinentes, para mandarlas al UAV inspector. El UAV inspector va a tener incorporado únicamente los datos y procedimientos necesarios para realizar la reconstrucción de viento. Con estas estimaciones en su mano, puede realizar modificaciones en su planificación, evitando zonas de riesgo, o variar la acción de control con previsión. Debido a que las operaciones de reconstrucción no requieren de alto coste computacional, se pueden realizar desde el propio UAV inspector, sin necesidad de hacer un paso previo con una estación de control en tierra. Con todo esto en mente, hay que abordar el primer problema para conseguir las estimaciones de viento, que es la construcción de un dataset. El dataset no es más que un conjunto de muchas simulaciones que contienen distintos casos de velocidades y orientaciones de viento. Hasta ahora esto no había sido necesario, ya que en los ejemplos del capítulo anterior 3.3.1, 3.3.2, el dataset ya estaba generado. 35 36 Capítulo 4. Estimación del Viento en Zonas de Inspección 4.1 Generación de Simulaciones Una vez se tiene clara la solución que se tiene entre manos, hay que llevarla a cabo. Primero, lo que se busca es representar el emplazamiento en un plano, es decir, en dos dimensiones, de tal forma que se tenga la planta del lugar. La zona de inspección va a constar de un edificio con dos depósitos, aunque la geometría es totalmente arbitraria, siendo posible llevar a cabo estructuras más complejas. El edificio va a ser rectangular de dimensiones 20m×6m , mientras que los depósitos van a ser circulares con un diámetro de 6m . Cabe destacar que no se busca una geometría complicada, ya que el objetivo es ilustrar que es posible obtener las predicciones de viento. Conocido el emplazamiento que se va a simular, hay que ver las distintas posibilidades de simulación. Lo que se busca es generar un dataset con muchas simulaciones donde haya variaciones de velocidad y orientación del viento. Estas variaciones no son totalmente aleatorias, debido a que el viento presenta ciertas actitudes regulares, en el sentido de que suele soplar en direcciones predominantes y sin llegar a superar umbrales de velocidad. Todo este tratamiento del aire, requiere ciertos conocimientos de mecánica de fluidos. Se va a suponer que el flujo es incompresible, es decir, número de Mach menor a la unidad, y que el régimen a simular es permanente, el transitorio en este caso no interesa. A su vez, hay que tener en cuenta el número de Reynolds, que depende del fluido, aire en este caso, y de la geometría donde se esté simulando. Como se verá, el número de Reynolds juega un papel crucial en la elección del solver. Con todos estos condicionantes a la hora de simular, hay que añadirle los parámetros que van a cambiar de simulación a simulación, velocidad del viento y ángulo con el que provenga. Lo que se busca es tener muchos casos, ya que conduce a un dataset rico en posibilidades. La solución que se obtendrá en cada simulación que compone el dataset, debe ser la velocidad del aire, en régimen permanente, en todo el emplazamiento que se está simulando. Por último, lo que queda es buscar distintos solvers de fluidos, que permitan llevar a cabo todas estas simulaciones. Se deberá hacer un estudio de las distintas posibilidades que pueden cumplir los requisitos para generar el dataset deseado. 4.1.1 Elección del Solver La elección del solver fluido es una parte crucial en el proceso de generar el dataset. Lo que se busca es un solver que sea capaz de mostrar buenos detalles en las simulaciones, pero que no requiera de un tiempo de ejecución excesivo. A su vez, hay que tener en cuenta que se va a trabajar en régimen turbulento, debido a que el aire por sus características inherentes presenta baja viscosidad, lo que implica un número de Reynolds alto. El hecho de intentar simular regímenes turbulentos conducirá a un solver más potente. Para empezar por algún solver, se pensó en que sería buena idea utilizar el solver empleado para generar el dataset del ejemplo del cilindro del capítulo anterior 3.3.3. En esta ocasión se entra más en detalle en que consiste este solver y las limitaciones que presenta. Estas limitaciones conducen a buscar un solver más fiable y potente. Mirando las distintas opciones, la que más se ajustaba a lo que se estaba buscando es Ansys Fluent. Immersed Boundary Method (IBM) El primer solver de fluidos, lo que implementa, es el método de la frontera inmersa, o Immersed Boundary Method (IBM) en inglés, [ 44 ], [ 46 ]. En particular, el método que se utiliza es el "método rápido", [ 45 ]. Este método lo que hace es resolver las ecuaciones de Navier-Stokes bidimensionales incompresibles alrededor de geometrías complejas. La principal característica que este tiene, es la capacidad para manejar geometrías complejas de manera eficiente, sin la necesidad de generar una malla de elementos finitos que se ajuste a la forma de los objetos sólidos, como puede hacer Ansys Fluent. A continuación, se procede a realizar una pequeña demostración de como sería simular el campo fluido haciendo uso de IBM. Una pequeña cuestión a tener en cuenta, es que estas simulaciones se hacen desde la ventana de comandos desde un sistema operativo LINUX, debido a la facilidad que este tiene para manejar comandos y añadirlos fácilmente. Las opciones de lanzamiento que permite aparecen en la figura 4.1. Verdaderamente, solo se hace uso de unas pocas variaciones. Por ejemplo, para realizar todo el dataset del ejemplo del cilindro 3.3.3, ejecutamos el siguiente comando: ../build/ibpm −nx 450 −ny 200 −ngrid 4 −length 9 −nsteps 4000 −xoffset −1−geom cylinder.geom −alpha 0 −Re 100 . Se entiende que si no se ha puesto algún argumento en la línea de comandos, se queda el valor que tenga por defecto. Por tanto, con estos parámetros se puede definir el dominio fluido. En caso de que se quiera implementar el mapa del 4.1 Generación de Simulaciones 37 Figura 4.1 Opciones de simulación de IBPM. emplazamiento, lo que habría que modificar es el archivo "cylinder.geom". Es en este archivo, donde se define la geometría de los distintos elementos que se encuentran dentro del dominio fluido. Aunque las posibilidades geométricas son bastante amplias, abarcando desde una simple línea hasta un perfil de ala, en la siguiente ilustración, figura 4.2, se mostrarán cuatro círculos que podrían simular cuatro tanques o silos: Figura 4.2 Representación de la vorticidad a través de cuatro cilindros. Con todo esto, se puede ver que es posible unas buenas simulaciones del fluido con un tratamiento sencillo de mallas y geometrías. Sin embargo, este método comienza a fallar cuando se ponen valores acordes a la realidad. Estos valores vienen a ser Reynolds muy altos, donde el fluido es el aire. 38 Capítulo 4. Estimación del Viento en Zonas de Inspección Para hacer posible un estudio lo más cercano a las condiciones reales, se había pensado que el fluido fuera aire y la velocidad del mismo no fuera superior a 10 m/s . Valores mayores que este, pondrían en peligro la seguridad del dron. Teniendo este valor de la velocidad, conocido el fluido que es aire y que la longitud característica puede considerarse el diámetro de un depósito, se estima que el número de Reynolds pertenecerá a un régimen turbulento muy alto, con un orden de 106. Bajo estas condiciones, se intenta probar el funcionamiento del método IBPM. Las simulaciones resultan erróneas, ya que en cada iteración para llegar a la solución, el error aumenta drásticamente. Al principio, no se conocía por qué el error aumentaba exponencialmente, pero tras hacer unos cálculos, se concluyó que se debe al mallado con el que se estaba trabajando. El objetivo, por tanto, era disminuir el mallado espacial, con los parámetros nx,ny yngrid, pero también disminuir el mallado temporal, parámetro dt. Sin embargo, tras muchos intentos y distintas pruebas que en el apartado teórico tenían sentido, no sé llegó a ninguna solución válida. Por tanto, se decidió a buscar otras posibilidades de solver. La clave de todo esto resulta en el valor tan grande de número de Reynolds con el que se pretende trabajar. Este solver funciona eficazmente y de manera fiable cuando los Reynolds son bajos, régimen laminar o régimen turbulento con valores bajos, valores que como máximo llegan a los 5000 de número de Reynolds. Es en ese momento cuando se plantea la posibilidad de utilizar Ansys Fluent. En un principio se descartó debido a que no se quería meter demasiada complejidad a la resolución de este problema. A su vez, este trabajo no tenía como fin hacer simulaciones de fluidos tan precisas como puede ofertar este programa. A pesar de todo esto, si se quería que las simulaciones fuesen de calidad, y más importante, que fuesen realizables, había que utilizar Ansys Fluent. Ansys Fluent Ansys Fluent es un software de simulación numérica de dinámica fluidos computacional. Es uno de los programas más utilizados en el campo de la investigación, debido a la alta calidad que presentan sus simulaciones. Es importante notar por qué se decidió a escoger este solver, y es principalmente porque se sabía de antemano que número de Reynolds tan altos, del orden de 106 , iban a ser tratados con solvencia. La simulación al completo se puede dividir en tres distintas partes: • Geometría. En esta primera parte lo que se hace es el diseño del emplazamiento como tal. La interfaz de Ansys a la hora de crear la geometría es bastante sencilla, lo que facilita poder implementar cualquier emplazamiento. En este caso se ha optado por un rectángulo y dos círculos, aunque podría haber sido algo más complejo. Cabe destacar que el dominio fluido ha de ser lo suficientemente grande comparado con la geometría en sí, para evitar que existan errores como flujos de recirculación. Aquí, figura 4.3, se puede ver cómo quedaría. • Mallado. El mallado es un mundo dentro de Ansys, por tanto, es importante saber qué es lo que buscas. En este caso es un ejemplo muy simple comparado con las posibilidades que puede llegar a ofrecer este mallado. Se toma de base la geometría realizada anteriormente y se divide el campo fluido en elementos finitos, generando el mallado. En esta ocasión se han elegido cuadriláteros de un tamaño máximo de siete metros. Lo más importante que se ha hecho en la malla es que cerca de las geometrías el mallado sea más fino, ya que será aquí donde surjan los grandes cambios, y no en zonas alejadas de la geometría. Con todo esto se obtiene el mallado de la figura 4.4. Por último, se puede determinar la calidad del mallado. Para ver esto, se puede usar el Orthogonal Quaility o el Skewness. Por un lado, un valor alto del Orthogonal Quaility significa que es una malla de calidad, es este caso la malla tiene un valor de media de 0.9872 . Por otro lado, el Skewness tiene que tener un valor bajo para conseguir una buena malla. En este caso, el valor del Skewness es 0.0567. • Configuración. Es el último paso antes de lanzar las simulaciones. Las simulaciones que se van a realizar van a ser en el permanente, ya que no interesa el transitorio en este caso, y el tipo de solver es basado en presión, debido a que no se va a superar el número de Mach. El modelo empleado de viscosidad es el modelo k-epsilon, aunque podría haberse escogido el modelo Spalart-Allmaras o el k-omega. El modelo k-epsilon se basa en la idea de descomponer el flujo turbulento en dos componentes principales: la energía cinética turbulenta (k) y la disipación de energía(epsilon). Estas dos variables se calculan para representar los efectos de la turbulencia. La elección de este modelo frente al Spalart-Allmaras es debido a que da soluciones más precisas a pesar de que sea computacionalmente un poco más costoso. El modelo k-omega da muy buenos resultados, en concreto en zonas de recirculación y en zonas críticas, sin embargo, es altamente costoso. Se concluyó que perder estos detalles de recirculación, no iban a ser relevantes en el resultado final. A su vez, es importante definir bien las condiciones de contorno. Se supone que el aire va a entrar con 4.1 Generación de Simulaciones 39 un ángulo máximo de ±20 ºy velocidades que van desde los 10 m/s hasta los 2m/s . Por tanto, como mínimo deberá de haber una entrada, que será del tipo velocity-inlet. A pesar de esto, se ha visto que el máximo de grados resulta un poco agresivo para las simulaciones, por lo que se van a poner dos entradas. Por otro lado, las salidas, que serán dos, son del tipo pressure-outlet. Para conseguir el dataset completo, el ángulo comenzará en 20 º, disminuyendo de dos en dos, hasta alcanzar los −20 º. Para cada ángulo habrá que examinar las nueve velocidades distintas, comenzando en 10 m/s y finalizando en 2m/s. Con todo esto, lo que se genera es un dataset de 189 casos distintos. Por último, dentro del solver basado en presión, el acoplamiento de presión-velocidad es resuelto mediante el algoritmo "SIMPLE". Los métodos de solución de discretización espacial usados serán todos de segundo orden, debido a la mayor precisión que estos presentan. El número de simulaciones es lo suficientemente alto para que los valores residuales sean bajos. Con todas estas consideraciones, se obtienen resultados como los de la figura 4.5. Figura 4.3 Representación de la geometría del emplazamiento. Figura 4.4 Representación del mallado. 40 Capítulo 4. Estimación del Viento en Zonas de Inspección Figura 4.5 Simulación de la velocidad absoluta del viento para una entrada de 10 m/s con inclinación nula. 4.2 Ubicación Óptima de Sensores y Reconstrucción En este último apartado, se verán los resultados de aplicar el algoritmo de ubicación óptima de sensores para el conjunto de datos generados con el solver explicado anteriormente. A su vez, se explicará detenidamente qué se ha hecho para obtener los resultados. 4.2.1 Generación de Journals Para empezar, hay que generar el set de datos. Como se ha explicado anteriormente, se obtendrán 189 simulaciones donde variará el ángulo de inclinación, dentro del rango ±20 º, y la velocidad de entrada, entre 10 m/s y 2m/s . Para un mejor entendimiento de esto último se realiza la ilustración de la figura 4.6. Estas simulaciones son obtenidas de manera automática, sin tener que ir lanzando cada simulación a mano. Esto es posible gracias a los journarls. Un journal es un archivo de secuencias de comandos y reproduce acciones realizadas por el ususario en el software de Ansys. Van a ser necesarios dos tipos de journals para generar todas las simulaciones. El primer journal, lo que hace es leer la geometría, el mallado y lanza la simulación en paralelo, cambiando en cada caso la velocidad de entrada en las componentes x e y . Como puede verse en la figura 4.6, tanto la parte superior funcionan como inlet youtlet. Esto es debido a que las simulaciones van a dividirse en dos bloques: el primero de ellos, recoge el intervalo de grados de 0 ºa 20 º, utilizando la parte superior como inlet y la inferior como outlet, mientras que el segundo bloque recoge el intervalo de grados de 0 ºa −20 º, utilizando la parte superior como outlet y la inferior como outlet. El motivo por el que se hace esto es causado por el alto nivel de inclinación que representan los ±20 º. Las simulaciones son lanzadas en paralelo, intentando que el máximo número de procesadores trabaje en la simulación. Una vez se lanza la simulación, este journal termina por guardar los archivos de simulación, de extensiones: .cas y.dat. El segundo journal realiza operaciones sencillas pero esenciales. Lo que hace es leer los archivos de simulación generados previamente por el otro journal,en serie, y guardar los resultados de la velocidad absoluta del aire en formato ASCII. El motivo por el que se crea un journal adicional, es debido a que si se guarda el formato ASCII en el mismo archivo donde se lanza la simulación en paralelo, los datos del archivo ASCII se encuentran corruptos. Esto es por problema que presenta Ansys al no permitir lanzar las 4.2 Ubicación Óptima de Sensores y Reconstrucción 41 simulaciones en paralelo y guardar los archivos en formato ASCII. Dentro de un archivo ASCII se encuentran las coordenadas x e y , junto a las magnitudes de la velocidad absoluta y la velocidad proyectada en las coordenadas x e y . La manera de almacenar los datos es tal que cada columna representa una simulación completa y en cada fila viene la información del fluido. La matriz de velocidades absolutas que se obtendrá es de 55084 ×189. Figura 4.6 Representación de las condiciones de contorno. 4.2.2 Aplicación de Algoritmos y Resultados Una vez que se tiene esa matriz, se procede a realizar las operaciones relacionadas con la reducción de dimensionalidad y la reconstrucción de datos: • Se comienza con la centralización de los datos en relación con la media. De esta manera, se elimina la influencia de la media y se enfoca en las variaciones respecto a ella. • Posteriormente, se escoge de manera aleatoria un determinado número de simulaciones para el entreno. En este caso se ha optado por escoger 150 simulaciones para construir el set de entrenamiento y 28 simulaciones para el set de validación. Es importante recalcar, que no se escogen las primeras 150 simulaciones, se escogen de manera aleatoria 150. Esto se hace para que el set de entrenamiento tenga en cuenta la mayoría de casos, propiciando la robustez. • Se realiza la descomposición de valores singulares de manera economizada al set de entrenamiento. Con esto se obtienen las matrices Ψ , S y V , de tal forma que: X=ΨSV −1 Se toma la matriz S , que contiene los valores singulares, para realizar el truncamiento óptimo de modos y sensores. Esto se representa gráficamente en la figura 4.7. El resultado de este truncamiento dará el número de hasta qué columna se considera la más representativa y enérgica de la matriz Ψ , pudiendo despreciarse el resto de las siguientes columnas. • Una vez que se sabe que el número de sensores mínimo para llevar a cabo la reconstrucción es 62, hay que aplicar los distintos algoritmos para encontrar la ubicación óptima. Todo esto viene explicado con más detenimiento en la sección del capítulo 3.2. En esta ocasión solo se utilizarán sensores aleatorios y sensores QR. Para ver en que cantidad se están cometiendo fallos, se realiza un estudio del error de reconstrucción. Este error de reconstrucción se calcula como el error relativo, donde se toma como medida verdadera una simulación del conjunto de validación. En el gráfico de la figura 4.8, hay que tener en cuenta varias cosas: la altura del rectángulo representa el rango del intercuartil; la mediana es la línea horizontal que se dibuja dentro del rectángulo; los valores atípicos u outliers, se representan mediante puntos individuales fuera del rango de los intervalos de confianza. Se puede extraer de la gráfica que conforme se añaden más modos se puede caer en sobreajuste y en 42 Capítulo 4. Estimación del Viento en Zonas de Inspección Figura 4.7 Representación de valores singulares normalizados. Se puede ver que se consigue el truncamiento óptimo con un valor de 62 modos. Esto implica que en el caso de que se tenga el mismo número de sensores que de modos, se colocarán 62 sensores. La imagen del modo 62 se obtiene al representar Ψ(62). una mayor incertidumbre, con lo que no conviene instalar un gran número de sensores. A su vez, la estrategia de sensores aleatorios no resulta favorable, debido al alto error relativo que se comete y la existencia de valores atípicos tan alejados de la mediana. Por último, se observa que el error cometido por los sensores QR es en torno al 10% , valor que se considera aceptable para la actividad que se tiene entre manos. 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 r 10-4 10-3 10-2 10-1 100 101 102 103 104 Error relativo Errores de Reconstrucción POD Sensores QR Sensores Aleatorios Figura 4.8 Errores cometidos cuando se lleva a cabo la reconstrucción. Para llevar a cabo la reconstrucción se utilizan los métodos de proyección sobre los r modos(base POD), sensores aleatorios y sensores que se obtienen con la factorización QR con pivotaje de columnas, subsección 3.2.1. 4.2 Ubicación Óptima de Sensores y Reconstrucción 43 • Finalmente, falta ver como queda esa reconstrucción y la ubicación de los sensores óptimos. Para ello lo que se hace es tomar una simulación del conjunto de validación y compararla con las reconstrucciones. Esto se hace en la ilustración de la figura 4.9. Por un lado, se puede ver que los sensores aleatorios se ubican por todo el dominio fluido para desencadenar en una reconstrucción pésima. Tan mala de reconstrucción, que existen zonas de velocidad negativa, cosa que no tiene sentido. Además, se puede ver que tiende al sobreajuste con solo unos pocos de sensores, como sucedió en la subsección de ejemplos 3.3.2. Por otro lado, se tiene que los sensores QR se disponen detrás de la estructura del emplazamiento. A su vez, para captar todos los detalles, en la parte inferior y superior derecha, se tienen unos pocos sensores que sirven para captar aquellos modos en los que en ángulo es extremo. La reconstrucción llevada a cabo por los sensores QR, es bastante acertada, en el sentido de que guarda perfectamente la forma que tiene el fluido y las velocidades existentes, a pesar de que deja de lado ciertos detalles. Hay que tener en mente que se han escogido 62 sensores de los 55084 posibles. Figura 4.9 Reconstrucción de una simulación aleatoria mediante la colocación de 62 sensores aleatorios y QR. 50 Bibliografía [48] K. Visalini, B. Subathra, S. Srinivasan, G. Palmieri, K. Bekiroglu, and S. Thiyaku, “Sensor placement algorithm with range constraints for precision agriculture,” IEEE Aerospace and Electronic Systems Magazine, vol. 34, no. 6, pp. 4–15, 2019. [49] Y. Guo, R. Lin, and M. D. Sacchi, “Optimal seismic sensor placement based on reinforcement learning approach: An example of obn acquisition design,” IEEE Transactions on Geoscience and Remote Sensing, vol. 61, pp. 1–12, 2023.