Algoritmo para la detección de formas aplicable a la estimación solar

Revista Iberoamericana de Autom´ atica e Inform´ atica Industrial 00 (2021) 1–11 www.revista-riai.org Algoritmo para la detecci´ on de formas aplicable a la estimaci´ on solar Aguilar-L´ opez, J.M.a,∗, Garc´ ıa, R.A.a, Camacho, E.F.a aDepartamento de Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Espa˜na. To cite this article: Aguilar-L´ opez, J.M., Garc´ ıa, R.A., Camacho, E.F. 2020. Two layer algorithm for spatial solar radiation estimation. Revista Iberoamericana de Autom´ atica e Inform´ atica Industrial 00, 1-5. https://doi.org/10.4995/riai.2020.7133 Resumen En este art´ ıculo se presenta un algoritmo h´ ıbrido bio-inspirado para la detecci´ on de formas aplicado a la estimaci´ on solar en plantas solares. Se tiene como objetivo localizar y caracterizar la forma de una nube sobre una planta solar bas´ andose en medidas de niveles bajos de la irradiancia con una peque˜ na flota de veh´ ıculos a´ ereos no tripulados (UAVs en ingl´ es) equipados con sensores capaces de medir la irradiancia directa normal. El algoritmo h´ ıbrido propuesto se inspira y adapta las ideas del algoritmo de optimizaci´ on de colonia de hormigas (ant colony optimization, ACO) y tambi´ en usa un algoritmo est´ andar de cobertura de ´ area, separ´ andose el campo de la planta solar en dos mallados, uno para cada capa del algoritmo, para encontrar el ´ area afectada por la nube. Cuando un UAV localiza la zona de baja irradiancia, los otros van a ayudarle. Dicho equipo delimita el borde de la nube usando conceptos de t´ ecnicas de procesamiento de im´ agenes. Finalmente, se prueba el algoritmo propuesto mediante simulaciones. Palabras clave: Estimaci´ on, Robots m´ oviles, Algoritmo de dos capas Shape detection algorithm applicable to solar estimation Abstract This paper presents a bio-inspired hybrid algorithm for shape detection applicable to solar estimation in solar power plants. The objective is to locate and characterise the shape of a cloud over a solar power plant based on low level irradiance measurement with a small fleet of Unmanned Aerial Vehicles (UAVs) equipped with direct normal irradiance sensors. The hybrid algorithm takes inspiration and adapts ideas of the ant colony optimisation algorithm (ACO) and also uses a standard cover area algorithm, separating the field into two grids, one for each layer of the algorithm, to find the area affected by the cloud. Once the low irradiance zone is located by one of the UAVs, the others go to help it. This team delimits the cloud border using concepts of an image processing technique. Finally, the algorithm is tested by simulations. Keywords: Estimation, Mobile robots, Two layers algorithm 1. Introducci´ on En los ´ ultimos a˜ nos, el avance de los veh´ ıculos a´ ereos no tripulados (unmmaned aerial vehicle, UAV) o drones, ha incentivado su uso en diferentes aplicaciones, tales como la agricultura (Rokhmana, 2015), las operaciones de salvamento (Silvagni et al., 2017), la prevenci´ on de desastres naturales (Yfantis, 2019), o la inspecci´ on de grandes ´ areas con problemas de accesibilidad por tierra (Cesetti et al., 2010). Las principales ventajes de este tipo de veh´ ıculos son sus reducidas dimensiones, maniobrabilidad, velocidad y capacidad de usar una gran variedad de dispositivos adaptados a la aplicaci´ on requerida, como c´ amaras RGB, t´ ermicas o hiperespectrales (Aasen et al., 2015), sensores de temperatura (Sheng et al., 2010), etc. Esta versatilidad los hace adecuados para un gran n´ umero de aplicaciones rob´ oticas. En rob´ otica, una de las l´ ıneas de investigaci´ on m´ as estudiadas es la conocida como Problema de Cobertura de ´ Areas o Problema de Cobertura de Camino (Galceran and Carreras, 2013). ∗Autor para correspondencia: [email protected] Attribution-NonCommercial-NoDerivatives 4,0 International (CC BY-NC-ND 4,0) Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 2 El objetivo es observar un ´ area de inter´ es tan r´ apido como sea posible, ya sea para tareas de vigilancia o para cualquier otro fin. Este problema es una variaci´ on del Problema del Viajante (Travelling Salesman Problem) (Current and Schilling, 1989). Existen numerosas aproximaciones para este problema, de las cuales algunas se basan en el tipo de descomposici´ on del ´ area, como por ejemplo, descomposiciones cl´ asicas en celdas como la de Morse (Acar et al., 2002) o la boustrofed´ on (Choset and Pignon, 1998); otras se pueden calcular fuera de l´ ınea (Oksanen and Visala, 2009) o en l´ ınea (Acar and Choset, 2000); y otras tienen en cuenta restricciones energ´ eticas (Choi et al., 2020). Para mejorar su eficiencia, normalmente los algoritmos para UAVs se desarrollan para sistemas multi-UAVs, es decir, para un equipo de UAVs. Este enfoque necesita algoritmos que controlen individualmente cada uno de ellos adem´ as de su formaci´ on conjunta (Anderson et al., 2008) (Garc´ ıa et al., 2020). La Naturaleza ha sido siempre una fuente de inspiraci´ on para la ciencia (Bar-Cohen, 2006), y esta clase de sistemas multiagentes se perfila para el empleo de algoritmos bio-inspirados como la Optimizaci´ on de Colonia de Hormigas (Ant Colony Optimisation ACO). Dicho algoritmo lo present´ o en 1991 Dorigo (Dorigo, 1991) como soluci´ on al Problema del Viajante (Johnson and McGeoch, 1997). Es un algoritmo bio-inspirado donde se crean hormigas virtuales que se mueven entre las distintas ciudades, representadas mediante nodos. Las hormigas dejan feromonas en los arcos que conectan los nodos. Al igual que las hormigas reales, las virtuales eligen su pr´ oximo nodo de destino bas´ andose en un m´ etodo probabil´ ıstico seg´ un la cantidad de feromonas que se ha dejado en el arco que conecta los nodos y la distancia entre ellos. Finalmente, tras cada iteraci´ on del algoritmo, una cantidad de las feromonas se evapora en todos los arcos. El algoritmo que se desarrolla en este trabajo ha sido probado en un ´ area de dimensiones similares a las que tendr´ ıa una t´ ıpica planta solar cilindropar´ abolica. Estas plantas ocupan grandes extensiones de terreno, normalmente 100 ha para plantas de 50 MW (S´ anchez et al., 2019). La eficiencia de una planta de este tipo puede verse tremendamente afectada por las sombras que causan las nubes (Ashley et al., 2017). En (Camacho et al., 2012; Camacho and Gallego, 2013) se muestra que es necesario saber la irradiancia directa normal (IDN o DNI en ingl´ es) para conseguir un control eficiente de la planta solar. Debido a que las plantas solares se extienden en grandes ´ areas (hasta 800 ha), algunas partes del campo solar pueden estar recibiendo mucha menos IDN por culpa del efecto de las nubes. Conocer la IDN en diferentes partes de la planta permite mejorar mucho el control de la misma (Xu et al., 2015). Esta IDN ha sido objeto de estudio en diferentes trabajos para mejorar la eficiencia de una planta, como (Nouri et al., 2018), o (Kuhn et al., 2017), donde se usan c´ amaras que apuntan al suelo para estimar la IDN. La principal desventaja de este tipo de sistemas es la necesidad de construir torres para la colocaci´ on de las c´ amaras y que, al ser una medida indirecta de la IDN obtenida a trav´ es de la reflectividad del suelo, se precisa de un modelo que depende de factores como la temperatura, humedad, la orograf´ ıa, etc., algunos de los cuales pueden variar a lo largo del a˜ no, lo que obligar´ ıa a reajustar el modelo. En (S´ anchez et al., 2018) se demuestra que con un control de los lazos adecuado mediante v´ alvulas se pueden obtener beneficios en la producci´ on de energ´ ıa, que se traducen en beneficios econ´ omicos. En el presente trabajo, para llevar a cabo la tarea de medir la IDN en varios puntos de la planta, se propone el uso de una flota de UAVs, equipado cada uno con un sensor liviano, de bajo consumo capaz de medir la IDN (Technologies, 2019). Este art´ ıculo presenta un novedoso algoritmo h´ ıbrido para la cobertura de ´ areas que combina un algoritmo inspirado en las ideas del algoritmo ACO para el problema de planificaci´ on de caminos y el algoritmo de boustrofed´ on para la cobertura del ´ area. El objetivo del algoritmo es detectar zonas con baja irradiancia por efecto de las nubes, y usar las medidas de la IDN en el campo para definir la forma de la sombra de la nube y para la estimaci´ on de la irradiancia espacial, en contraposici´ on a lo presentado en otros trabajos recogidos en (Law et al., 2014). La principal contribuci´ on de este trabajo respecto a otros, como por ejemplo (Xiong et al., 2019), es su car´ acter h´ ıbrido y su aplicaci´ on a la detecci´ on de ´ areas con poca irradiancia en plantas solares. El resto del art´ ıculo se organiza como se indica a continuaci´ on. En la Secci´ on 2 se presentan la notaci´ on y los conceptos preliminares relacionados con el algoritmo ACO y la descomposici´ on del terreno. En la Secci´ on 3 se describe el problema que se ha considerado. En la Secci´ on 4 se muestra el algoritmo h´ ıbrido propuesto, que se pone a prueba por simulaci´ on en la Secci´ on 5. Finalmente, las conclusiones se detallan en la Secci´ on 6, donde tambi´ en se remarcan ciertas consideraciones. 2. Notaci´ on y preliminares 2.1. Algoritmo ACO original Como se ha mencionado en la secci´ on anterior, este algoritmo fue inicialmente propuesto por Dorigo (Dorigo, 1991). A continuaci´ on, se describe el algoritmo particularizado para la resolucion del Problema del Viajante. Un conjunto de hormigas B={h1,h2,...,hi, . . . , hH}tienen que conectar el hormiguero, que es el nodo de partida, al lugar con la comida, que es el nodo de destino. Ambos nodos pertenecen al conjunto de nodos V={n1,n2,...,ni, . . . , nN}. Las hormigas completan su misi´ on viajando a trav´ es de los arcos que conectan los nodos, depositando feromonas cuando dejan dichos arcos. Las hormigas eligen de forma probabil´ ıstica cu´ al ser´ a el siguiente nodo que visitan, seg´ un la cantidad de feromonas que haya en el arco de uni´ on entre su nodo y los nodos conectados a ´ el y seg´ un la distancia entre los nodos. Sea di j la distancia entre el nodo iy el nodo j, se define vi j la visibilidad del nodo idesde el nodo jpara el caso del Problema del Viajante como: vi j =1 di j (1) La probabilidad de que una hormiga hse mueva desde el nodo ial nodo jsiendo τi j la cantidad de feromonas en el arco que conecta el nodo icon el nodo jes: ph i j =[τi j(k)]α[vi j]β PgJh i[τig(k)]α[vig]β,∀jJh i, α ≥0, β ≥0,(2) donde Jh ison los nodos no visitados por la hormiga hdesde el nodo i. Los par´ ametros αyβajustan el peso relativo de la Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 3 heur´ ıstica en el c´ alculo de la probabilidad. El caso de α=0 es el t´ ıpico algoritmo codicioso o voraz, donde los nodos m´ as cercanos son los m´ as probables de visitar (Zhang et al., 2000), mientras que en el caso β=0 solo se considera el peso de las feromonas. La evaporaci´ on de las feromonas en cada iteraci´ on se calcula con Eq. (3): τi j(k+1) =(1 −ρ)τi j(k)+ H X h=1 ∆τh i j(k),(3) siendo ρla variable que ajusta el grado de evaporaci´ on y ∆τh i j(k) la cantidad de feromonas que deposita la hormiga hen el arco i j en el instante k. Finalmente, el algoritmo ACO usa una funci´ on objetivo que hay que maximizar o minimizar seg´ un el tipo de problema. Particularmente, en el Problema del Viajante la funci´ on es: m´ ın N X i=0 N X j,i,j=0 ci j xi j.(4) siendo ci j el coste de viajar por el arco que conecta los nodos iyj, y xi j una variable binaria que indica si se ha viajado por el arco i j. En resumen, el algoritmo ACO consta de los siguientes pasos: 1. Las hormigas empiezan en nodos aleatorios. 2. Cada hormiga hcrea un camino soluci´ on Phentre su nodo de partida y el de destino visitando los otros nodos y eligiendo el arco que los conecta de forma probabil´ ıstica, seg´ un Eq. (2). 3. Se eval´ uan todos los caminos soluci´ on Phcon la funci´ on objetivo Eq. (4) y cada hormiga hdeja una cantidad de feromonas en los arcos que ha recorrido proporcional a la calidad de su soluci´ on. 4. El mejor camino soluci´ on encontrado en esta iteraci´ on Ph se compara con el mejor camino soluci´ on guardado, y si lo supera, lo sustituye. 5. Se evapora la cantidad de feromonas de cada arco i j seg´ un Eq. (3). 6. Se empieza de nuevo desde el paso 1. Para decidir si el ´ ultimo camino soluci´ on almacenado es una soluci´ on aceptable, se establece una condici´ on de fin que se eval´ ua despu´ es de cada iteraci´ on, ya que no hay otra forma de saber si la soluci´ on encontrada es el m´ aximo o m´ ınimo absoluto. 2.2. Descomposici´on del campo Se divide el terreno en dos capas, una para el algoritmo inspirado en ACO propuesto y otra para el algoritmo de cobertura de ´ area. Usando la terminolog´ ıa de jerarqu´ ıas, el mallado para el primer algoritmo se denomina malla superior y el mallado para el segundo malla inferior. La Figura 1 muestra ambas mallas: la malla superior se representa en azul y la inferior en rojo. La malla superior del terreno se divide en un n´ umero entero de celdas del mismo tama˜ no, 10×10 en este caso, siendo cada una de 100 ×100 metros. Esta descomposici´ on se elige tras comprobar que una divisi´ on mayor en el mallado superior no aprovechar´ ıa las ventajas del algoritmo de cobertura de ´ area, y uno m´ as peque˜ no no aprovechar´ ıa las ventajas del algoritmo inspirado en ACO, incrementando el tiempo de b´ usqueda final. De igual forma, la malla inferior tambi´ en se divide en un n´ umero entero de celdas. El detalle de esta descomposici´ on puede observarse en la Figura 1, donde se observa una divisi´ on de 5 ×5 celdas, siendo cada una de 20 ×20 metros. La descomposici´ on de la malla inferior mostrada es diferente a la que se usa en el algoritmo para que se pueda distinguir en la figura, ya que en las simulaciones realizadas cada celda es de 1×1 metros, dando, por tanto, una malla de 1000×1000 celdas. Figura 1: Esquema del mallado de las dos capas usadas en el algoritmo propuesto. 2.3. Sistema multi-UAVs El modelo de comunicaci´ on del sistema multi-UAVs utilizado es el de comunicaci´ on centralizada: los UAVs no se comunican directamente entre s´ ı, sino que env´ ıan sus medidas y reciben ´ ordenes del ordenador de la estaci´ on base, que es el que procesa los datos y calcula las nuevas rutas. La principal desventaja de este tipo de comunicaci´ on es la necesidad de un canal seguro de comunicaci´ on entre los sistemas que lo componen. La aplicaci´ on de esta propuesta est´ a orientada a plantas solares que generalmente tienen una extensi´ on de 1 km2, permitiendo las redes Wi-Fi actuales garantizar esta comunicaci´ on en estas distancias y a una velocidad adecuada. Otro de los puntos a tratar cuando se trabaja en un sistema multi-robot como el propuesto es c´ omo evitar la colisi´ on entre los distintos robots del sistema. La soluci´ on general pasa por el uso de sensores junto a una buena comunicaci´ on entre los robots o algoritmos de generaci´ on de trayectorias que tengan esto en cuenta. En el problema presentado, donde se usan UAVs, se pueden evitar f´ acilmente las colisiones estableciendo distintas alturas de seguridad para cada uno de los UAVs. Para poder medir la irradiancia solar, los UAVs disponen, como ya se mencion´ o en la Secci´ on 1, de un sensor ligero y de bajo consumo, concretamente de 10 gramos de peso y un campo de visi´ on de hasta 120º. Aunque este trabajo es te´ orico y se prueba en simulaci´ on, cuando en futuros trabajos se pruebe experimentalmente, se garantizar´ a que la medida no se vea afectada por una mala orientaci´ on del campo de visi´ on del sensor utilizando un Gimbal que lo oriente en los ´ angulos deseados. La conjunci´ on del gran campo de visi´ on del sensor y del Gimbal permitir´ a obtener buenos resultados a´ un con peque˜ nos desplazamientos o vibraciones de los UAVs. Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 4 3. Planteamiento del problema En esta secci´ on se describe el problema propuesto. El objetivo del equipo de UAVs es localizar y caracterizar la forma de una nube Cpara mejorar la eficiencia de una planta solar situada en un ´ area delimitada A. El equipo est´ a formado por un conjunto de UAVs denotados como U={u1,u2, . . . , ui, . . . , uH}. Como se ver´ a luego, cada UAV requiere su propia hormiga virtual y, por tanto, el n´ umero de hormigas es igual al n´ umero de UAVs. El ´ area de inter´ es A es un rect´ angulo de wpor l. Las siguientes consideraciones se han tenido en cuenta para la soluci´ on propuesta, que es un primer acercamiento al problema planteado. Se aproxima la nube Cpor una elipse con semieje mayor a, semieje menor b, centro Oy´ angulo de rotaci´ on φ. V´ ease la Figura 2. Dado que el objetivo es detectar nubes de tama˜ no suficiente para que afecten a la eficiencia de la planta solar, se establece un semieje mayor y menor m´ ınimo de 20 metros. Se supone que la nube se mueve tan despacio que puede considerarse est´ atica. La reducci´ on de la irradiancia solar provocada por la nube Cse modela como una funci´ on lineal, v´ ease la Figura 3. Cse divide en tres subregiones con diferente factor de oclusi´ on, como se aprecia en la Figura 4. Estas subregiones tambi´ en son el´ ıpticas y mantienen la misma excentricidad que Cpero tienen diferente semieje mayor. En futuros trabajos muchas de estas restricciones se ver´ an relajadas para abarcar un espectro m´ as amplio de situaciones, como nubes no convexas o en movimiento, que ser´ an comentadas en la Secci´ on 6. Figura 2: Forma de la nube supuesta. 4. Soluci´ on propuesta El algoritmo h´ ıbrido propuesto se descompone en dos fases: 1. La b´ usqueda de la nube: el objetivo de esta fase es encontrar una nube en el ´ area A. Para ello, los UAVs miden continuamente la irradiancia solar, que decaer´ a bajo la influencia de la sombra de una nube. Esta primera fase se divide, a su vez, en dos etapas: a)Etapa del algoritmo inspirado en ACO: cada UAV decide cu´ al ser´ a la siguiente celda de la malla superior a la que se dirigir´ a (ver Subsecci´ on 4.1). b)Etapa de cobertura de la celda: cada UAV barre la celda de la malla superior escogida visitando parte de las celdas de la malla inferior contenidas en ella (ver Subsecci´ on 4.2). Si uno de los UAVs toma una medida por debajo de un determinado valor umbral, se considerar´ a que se ha encontrado una nube y concluir´ a esta fase; en caso contrario, se contin´ ua alternando entre estas dos etapas. Se notificar´ a este hecho a los otros UAVs para que se dirijan a esta posici´ on y se pase a la siguiente fase. En caso de que en el camino a dicha posici´ on detecten otro punto distinto de la nube, se quedar´ an en esa posici´ on para continuar con la siguiente fase desde ah´ ı. Este valor umbral no es un valor fijo, y se puede variar durante la ejecuci´ on del algoritmo. Es un valor de referencia que se escoge bas´ andose en la IDN del cielo despejado y el error de los sensores. El valor de la IDN se puede conocer tanto por datos meteorol´ ogicos de a˜ nos anteriores como a trav´ es de una medida de un sensor que apunte directamente al sol, pudi´ endose usar el mismo sensor que montan los UAVs u otro m´ as preciso. Figura 3: Funci´ on de irradiancia supuesta, midiendo Pro f undidad de la nube la distancia desde el borde de la nube hasta su centro, es decir, pro f undidad =0 es el borde de la nube. Figura 4: Densidad de la nube supuesta. Hay dos subregiones dentro de la nube: la primera es la mostrada en azul claro y la segunda en azul oscuro. El centro de la nube, por ser la zona m´ as profunda, tiene el valor m´ ınimo de irradiancia Imin y es por ello la zona m´ as oscura. El valor m´ aximo de irradiancia Imax est´ a representado por el color amarillo anaranjado del borde. 2. Determinar la forma de la nube: en esta fase los UAVs encuentran el borde de la nube y lo recorren. Primero calculan el gradiente de irradiancia solar. Despu´ es se mueven siguiendo el vector perpendicular a dicho gradiente, es decir, a lo largo de la isol´ ınea de irradiancia solar, corrigendo el vector de movimiento poco a poco cuando Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 5 detectan que se salen de la nube a trav´ es de las medidas de irradiancia solar tomadas. Con este movimiento a trav´ es de la isol´ ınea se consigue delimitar el contorno de la nube, y con ello, su forma. Esta fase se divide en tres etapas: a)Medida en espiral cuadr´atica de la irradiancia solar: al principio no se puede calcular el gradiente porque solo se ha medido un punto de la nube. Por ello, cada UAV describe una espiral cuadr´ atica (como la espiral de Ulam (Daus, 1932), v´ ease la Figura 7) para tomar suficientes medidas y poder calcular el gradiente (ver Subsecci´ on 4.3). b)C´alculo del gradiente de irradiancia solar: para calcular el gradiente se usan las medidas tomadas en la etapa anterior con el operador Sobel, una t´ ecnica de p´ ıxeles de procesamiento de im´ agenes y de visi´ on por computador (ver Subsecci´ on 4.4). c)Movimiento a lo largo de la isol´ınea de irradiancia solar: los UAVs se mueven siguiendo el vector perpendicular al gradiente calculado, manteniendo un sentido antihorario, aunque hubiera sido indiferente haberlo hecho en sentido horario. En el momento en el que detecten que est´ an fuera de la zona de influencia de la nube, girar´ an suavemente hacia la izquierda los grados necesarios para volver a encontrar el borde de la nube (ver Subsecci´ on 4.5). En las siguientes subsecciones se describen cada una de estas fases y etapas. 4.1. Algoritmo inspirado en ACO para cobertura de ´areas El problema abordado en este trabajo no usa el algoritmo ACO como una optimizaci´ on tal y como se ha visto en la Secci´ on 2, sino que se inspira en sus ideas para resolver la cobertura de ´ area. Para ello se han hecho algunas modificaciones a la adaptaci´ on del comportamiento de las hormigas presentadas en el algoritmo original. Cabe se˜ nalar que la descomposici´ on en celdas del mallado superior es el equivalente a los nodos del grafo asociado a un problema de optimizaci´ on resuelto con un algoritmo ACO. En la explicaci´ on del algoritmo propuesto se hablar´ a de celdas y nodos indistintamente. Para empezar, para una cobertura r´ apida del terreno, las hormigas deben ser repelidas por las feromonas, en lugar de ser atra´ ıdas por ellas, haci´ endose as´ ı m´ as probable que se visiten nodos donde haya menos feromonas. Por tanto, la probabilidad de escoger un nodo jviene dada por: pτ i j =1−τj Pτg ,(5) donde Pτges la cantidad de feromonas de todos los nodos. Otro cambio importante es el que se ha realizado en el criterio de la distancia. En el algoritmo ACO original, para el Problema del Viajante, se trata con el c´ alculo de la visibilidad Eq. (1): el nodo m´ as lejano tiene menos probabilidad de ser escogido. En el caso de la cobertura de ´ area, el objetivo es visitar cada celda al menos una vez, sin importar la distancia que haya que recorrer. Por tanto, en este caso, el criterio de la distancia no se aplica. Alternativamente, para cubrir el ´ area m´ axima posible, es una ventaja que las hormigas no est´ en cercanas entre s´ ı, por lo que el criterio de distancia se sustituye por un criterio de distancia m´ ınima con respecto a otras hormigas. Esto posibilita que las hormigas se repelan entre s´ ı y no visiten a la vez el mismo nodo. (a) C´ ırculo umbral R: se evita la hormiga h2mientras que la hormiga h3no es tenida en cuenta ya que se encuentra fuera del c´ ırculo R. Las flechas azules muestran algunos movimientos posibles para evitar a la hormiga h2. (b) A efectos de claridad en la representaci´ on solo se muestran algunas distancias d, aunque el algoritmo las eval´ ua todas. Primero se calculan las distancias entre la celda a evitar, la de h2, y las candidatas (i,j,k). A continuaci´ on, se le asigna a cada celda una probabilidad inversamente proporcional a su distancia drespecto a la celda a evitar. La celda rayada es la m´ as probable para ser escogida puesto que est´ a m´ as lejos de la celda de h2. Figura 5: Criterio de distancia con otras hormigas. El criterio de distancia con respecto a otras hormigas funciona de manera similar al criterio de las feromonas. Solo se le aplica a aquellas hormigas hjque est´ en por debajo de una cierta distancia umbral, es decir, las que se encuentran dentro del c´ ırculo de radio R. Sea una hormiga h1en el centro de un c´ ırculo de radio Rseparada de otras hormigas h2yh3una distancia rh2 yrh3respectivamente. En la Figura 5(a) se puede observar que si r<R, entonces la hormiga en cuesti´ on se evita, en el ejemplo h2, mientras que en otro caso se ignora, en el ejemplo h3. Una nueva probabilidad se asigna a cada celda, que ser´ a mayor cuanto m´ as lejos est´ e de la hormiga h2: pD i j =dj Pdl ,(6) siendo djla distancia entre la celda candidata, celda j, y la celda ocupada por la hormiga h2, como se muestra en la Figura Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 6 5(b). Adem´ as, dles la distancia entre la celda de h2y cualquier otra celda candidata por lo que Pdlel sumatorio de las distancias entre la celda ya ocupada por h2y todas las posibles. En caso de que hubiera m´ as de una hormiga a evitar, se calcular´ ıa independientemente cada pi j respecto de cada una de ellas y se aplicar´ ıa la ecuaci´ on (7): pD i j = pD1 i j pD2 i j PpD1 i j pD2 i j (7) Finalmente, la probabilidad de que la hormiga h1escoja una celda como nuevo destino se obtiene combinando el criterio de las feromonas y el criterio de distancia con otras hormigas, seg´ un la ecuaci´ on: ph1 i j = pτ i j pD i j Ppτ ig pD ig .(8) El algoritmo propuesto elige una celda del mallado superior como siguiente posici´ on de h1una vez calculadas todas las probabilidades de las celdas candidatas. 4.2. Algoritmo de cobertura de celdas Una vez escogida la nueva celda del mallado superior, el algoritmo de cobertura crea un camino compuesto de puntos (waypoints) para recorrer las celdas del mallado inferior contenidas en ella. Hay que destacar que no se recorren todas las celdas del mallado inferior, ya que solo se quieren detectar nubes mayores de un determinado tama˜ no m´ ınimo, par´ ametro que se tiene en cuenta para dise˜ nar el recorrido que barre las celdas. Se dise˜ naron dos tipos de caminos de cobertura: (a) Pasillos. (b) Diagonal. Figura 6: Caminos de cobertura usados en el algoritmo h´ ıbrido propuesto. 1. Pasillos: el UAV se mueve desde un lateral de la celda al otro, se desplaza entonces perpendicularmente una distancia concreta, y vuelve al lateral desde el que empez´ o. Esto se repite hasta que se cubra completamente la celda del mallado superior, trazando un camino que se asemeja a varios pasillos. Este tipo de movimiento se conoce como boustrofed´ on, explicado con m´ as detalle en la Subsecci´ on 5.2. El ancho del pasillo, o paso, se fija seg´ un el tama˜ no m´ ınimo de nube que se quiere detectar. Como en este caso el m´ ınimo de los semiejes es de 20 metros (ver Secci´ on 3), ese ser´ a el ancho del pasillo. Un ejemplo de este tipo de camino se muestra en la Figura 6(a). Este camino tiende a ser generalmente m´ as lento pero tambi´ en m´ as exhaustivo que la otra propuesta, el camino diagonal. A pesar de todo, la velocidad depende en gran medida del tama˜ no de la celda y el paso que se use, siendo m´ as r´ apido cuanto m´ as peque˜ na sea la celda o cuanto m´ as ancho sea el paso. 2. Diagonal: este camino solo recorre las diagonales de la celda y dos de los lados laterales, por lo que normalmente es m´ as corto que el camino de pasillos, ya que se visitan menos celdas del mallado inferior, aunque depende del ancho del pasillo escogido. No obstante, deja sin visitar una superficie de la celda considerable. Observaci´ on: tras probar ambos tipos de camino se ha descartado el camino tipo diagonal debido a la superficie sin explorar que deja este camino. Cuando se termina de explorar la celda del mallado superior, se vuelve al algoritmo inspirado en ACO para elegir una nueva celda y cubrirla con un nuevo camino. Este proceso se repite hasta que alg´ un UAV detecte un punto con baja irradiancia, es decir, un punto de una nube, y llame a los dem´ as para empezar la fase que definir´ a la forma de la nube. 4.3. Medida de la irradiaci´on solar mediante la espiral cuadr´atica Cuando se detecta la nube solo se conoce un punto de ´ esta, insuficiente para calcular el gradiente de la irradiaci´ on solar. Para obtener m´ as medidas se usa un camino en forma de espiral cuadr´ atica como la que se ve en la Figura 7. Cuanto mayor sea la espiral, mejor ser´ a el c´ alculo del gradiente, aunque tambi´ en har´ a m´ as lento el algoritmo pues al UAV le lleva m´ as tiempo recorrer el camino. Se establece as´ ı una relaci´ on de compromiso entre el la calidad del c´ alculo y la velocidad del algoritmo. Para determinar el n´ umero de giros utilizados, se han realizado 1000 simulaciones de medidas de espirales de hasta 5 giros con un 3 % de ruido, mostr´ andose los resultados en la Tabla 1, siendo el error la desviaci´ on en grados sexagecimales del vector gradiente obtenido respecto del vector gradiente real, que es el que se obtendr´ ıa si se conociera la irradiancia exacta o la funci´ on del factor de oclusi´ on de la nube. Tabla 1: Simulaciones de la espiral. N´ umero de giros Error (º) Tiempo (s) 1 43.68 1.05 2 19.15 3.2 3 9.61 6.12 4 6.17 10.31 5 4.13 15.47 A ra´ ız de los resultados, se fij´ o el n´ umero de giros en 3, dado que 10º es una desviaci´ on asumible. Este valor de 10º se determin´ o a trav´ es de pruebas de ensayo y error con el algoritmo. En cualquier caso, este n´ umero de giros podr´ ıa determinarse en ejecuci´ on en funci´ on de la distancia a la planta, dado que si est´ a m´ as lejos de la misma podr´ ıa emplear m´ as tiempo para obtener m´ as precisi´ on. Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 7 4.4. Gradiente de irradiancia solar El gradiente de irradiancia solar se calcula usando el operador Sobel, un algoritmo de p´ ıxeles que se utiliza en el procesamiento de im´ agenes en el campo de la visi´ on por computador, desarrollado por Irwin Sobel y Gary Feldman en 1968 (Sobel and Feldman, 1968). Es una diferenciaci´ on discreta usada principalmente para detectar bordes en las im´ agenes. Figura 7: Camino de espiral cuadr´ atica usado para tomar medidas de la nube. El algoritmo usa dos kernels 3x3, uno para la direcci´ on horizontal y otro para la vertical. Las ecuaciones Eq. (9a) y Eq. (9b) muestran los kernels que se usan en el algoritmo de Sobel. Posteriormente, Prewitt, Scharr o Roberts Cruz propusieron otros kernels alternativos para obtener resultados diferentes (Savant, 2014). En cualquier caso, en este trabajo solo se han utilizado los originales. KX=          −101 −202 −101           ,(9a) KY=          121 000 −1−2−1           .(9b) Como se ve en Eq. (9a) y Eq. (9b), la aproximaci´ on discreta de las derivadas en los ejes xeyse calcula convolucionando los kernels con el mapa de irradiaci´ on solar I. Este mapa se puede considerar como una imagen de tantos p´ ıxeles como celdas del mallado inferior se hayan visitado en el camino de espiral cuadr´ atica, siendo as´ ı la intensidad de la irradiaci´ on en cada celda la intensidad de color del p´ ıxel correspondiente en la imagen. Cada p´ ıxel en Ise puede etiquetar de dos maneras: Evaluado: uno de los UAVs ha visitado la celda del mallado inferior equivalente y ha medido su valor de irradiancia. Este valor es un float positivo de dos decimales. Estado desconocido: ning´ un UAV ha visitado la celda del mallado inferior correspondiente. Su valor es -1. Figura 8: Campo de gradientes y el vector gradiente medio resultante despu´ es de la exploraci´ on de la nube con el camino de espiral cuadr´ atica. El m´ odulo del vector gradiente medio est´ a amplificado a efectos de una representaci´ on clara. En general, antes de comenzar el barrido en espiral cuadr´ atica, los p´ ıxeles de Iest´ an marcados con −1, salvo por aquellas localizaciones que previamente hayan sido visitadas por un UAV durante las fases anteriores del algoritmo. Una vez que acaba el recorrido en espiral, todos los p´ ıxeles han sido evaluados y est´ an listos para que se les aplique el operador de Sobel. El gradiente de irradiancia solar se obtiene con Eq. (10a) y Eq. (10b). Ix=δI δx=KX∗I,Iy=δI δy=KY∗I,(10a) ∇I="Ix Iy#.(10b) Usando el operador Sobel sobre todo el mapa Ise obtiene un campo de gradientes. Calculando la media de todos los vectores gradiente del campo se obtiene el vector gradiente medio desde la posici´ on inicial en el borde de la nube, como se muestra en la Figura 8. Este vector se usa para determinar la direcci´ on de movimiento del UAV. 4.5. Movimiento siguiendo la isol´ınea Teniendo el gradiente resultante, la isol´ ınea de irradiancia solar est´ a en la direcci´ on perpendicular al mismo y es la que sigue el UAV en sentido antihorario para recorrer el borde de la nube. Durante este movimiento el vector de trayectoria se vuelve r´ apidamente impreciso ya que la isol´ ınea es el´ ıptica. Este error se corrige continuamente, ya que, cuando el UAV detecta que est´ a fuera de la nube, a trav´ es de una medici´ on de la irradiancia por encima del valor umbral, gira su trayectoria suavemente en sentido antihorario para volver a la nube. De esta manera la forma de la nube se delimita paso a paso. Cabe se˜ nalar que cuanto m´ as lento se muevan los UAVs menos tardar´ an en corregir la trayectoria cuando se salgan de la nube, pero tambi´ en les llevar´ a m´ as tiempo definir la forma completa de la misma. La Figura 9 muestra el borde que se obtiene de la nube y tambi´ en c´ omo los 3 UAVs modifican su trayectoria con el procedimiento explicado anteriormente. Se puede apreciar que las l´ ıneas son un poco onduladas debido a la correcci´ on descrita. Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 8 Figura 9: Forma de la nube definida por los 3 UAVs (azul, verde, morado) aplicando el algoritmo propuesto. 5. Simulaciones Para probar el m´ etodo h´ ıbrido propuesto se dise˜ naron una serie de simulaciones. ´ Estas fueron programadas y ejecutadas en el programa MATLAB versi´ on R2020a. Para una mejor representaci´ on gr´ afica de los resultados las figuras de esta secci´ on muestran un ´ area y una nube que son m´ as peque˜ nas que las usadas en las simulaciones, pero proporcionales a ´ estas. (a) Ejemplo de b´ usqueda de la nube usando el m´ etodo h´ ıbrido propuesto. (b) Ejemplo de b´ usqueda de la nube usando el movimiento boustrofed´ on. Figura 10: B´ usqueda de la nube usando los dos algoritmos. 5.1. Condiciones de las simulaciones Cada simulaci´ on es un escenario de b´ usqueda de la nube en un ´ area de 1000x1000 metros. La descomposici´ on de la superficie fue de 10x10 celdas para el mallado superior (algoritmo ACO), que es visible en las figuras, y de 1000x1000 celdas para el mallado inferior (movimiento boustrofed´ on), que no se muestra en las figuras por su peque˜ no tama˜ no. El equipo de UAVs est´ a formado por 3 veh´ ıculos, y cada uno de ellos se mueve a 8 m/s sin restricciones de movimiento ni de energ´ ıa. El valor de la velocidad se basa en la velocidad media de UAVs comerciales como (DJI Technology Inc., 2015). Sus posiciones iniciales, en metros, son: (0,333.33), (1000,666.66) y (0,1000). La nube se modela como una elipse de 100 metros de semieje mayor y 80 metros de semieje menor. En cada simulaci´ on, los par´ ametros que var´ ıan son la posici´ on y el ´ angulo de giro de la nube, que se generan aleatoriamente cada vez. Aguilar-L´opez, J.M. et al. /Revista Iberoamericana de Autom´atica e Inform´atica Industrial 00 (2021) 1–11 9 (a) Algoritmo h´ ıbrido (b) Algoritmo boustrofed´ on. Figura 11: Cuando se encuentra la nube se llama al resto de UAVs. Las l´ ıneas continuas representan el camino que se ha seguido antes de encontrar la nube, mientras las discontinuas muestran el camino hacia la nube tras la llamada. Adem´ as, cada escenario se usa para evaluar tanto el algoritmo h´ ıbrido como el algoritmo boustrofed´ on cl´ asico, ya que ´ este es uno de los m´ as usados en la literatura para abordar el problema de cobertura de ´ area, como se menciona en la Secci´ on 1. 5.2. Algoritmo boustrofed´on Los algoritmos de tipo boustrofed´ on son los que est´ an inspirados en un tipo de textos bidireccionales, esencialmente textos antiguos, en los que la direcci´ on de la escritura de cada l´ ınea es contraria a la direcci´ on de la l´ ınea anterior, es decir: la primera l´ ınea se escribe de izquierda a derecha, la segunda de derecha a izquierda y as´ ı sucesivamente. Este mismo procedimiento se aplica para generar caminos que cubran un ´ area, pues asegura que todas las posiciones del ´ area se visiten dentro de un cierto tiempo usando un camino bastante eficiente sin demasiados giros, como puede verse en la Figura 10(b). Algunos algoritmos inspirados por este tipo de movimiento se pueden ver en (Coombes et al., 2017), (Ntawumenyikizaba et al., 2012), (Jin and Tang, 2010) or (Avellar et al., 2015). El camino de pasillos del m´ etodo h´ ıbrido es un ejemplo de este tipo de algoritmos boustrofed´ on, como se ve en la Figura 10(a). Es importante subrayar que la implementaci´ on del algoritmo boustrofed´ on utilizada aqu´ ı para compararse con el algoritmo h´ ıbrido no visita todas las celdas del ´ area por la raz´ on explicada en la Subsecci´ on 4.2, esto es, solo las nubes de un determinado tama˜ no m´ ınimo o superior son de inter´ es. Por esta raz´ on, este algoritmo es, en realidad, como si se usara el camino de pasillos que se usa para cubrir una celda del mallado superior, pero considerando toda el ´ area como una ´ unica celda. El ancho de separaci´ on entre pasillos es el mismo que se usa en el algoritmo h´ ıbrido, es decir, 20 metros. Adem´ as, para evitar solapamientos, cada UAV tiene asignado un tercio del ´ area, empezando su exploraci´ on en una esquina, de forma que se optimiza su movimiento y b´ usqueda en el ´ area asignada. Como en el caso del algoritmo h´ ıbrido, los UAVs que aplican este algoritmo est´ an continuamente midiendo la irradiancia para detectar la nube. En la Figura 11 se ve c´ omo, gracias a esta continua medici´ on, al encontrar uno de los UAVs una medida por debajo del umbral los dem´ as son llamados para continuar con la fase de descripci´ on de la forma de la nube. Figura 12: Comparativa de los resultados de las simulaciones de ambos algoritmos. 5.3. Resultados Los resultados obtenidos tras ejecutar 1000 simulaciones se muestran en la Tabla 2. Cabe se˜ nalar que se compara la velocidad en la detecci´ on de la nube, y no su posterior caracterizaci´ on, dado que para esa fase ambos algoritmos emplean el mismo m´ etodo. El algoritmo h´ ıbrido propuesto obtiene mejores resultados que el boustrofed´ on en t´ erminos generales, aunque con margen de mejora: le lleva un 21 % menos de tiempo de media encontrar la nube. Es importante se˜ nalar la existencia de valores at´ ıpicos, como se puede ver en la Figura 12. Este tipo de valores empeoran el valor de la media, pero la mediana es bastante m´ as robusta a ellos, raz´ on por la que se recoge en la Tabla 2. Un