Simulación dinámica de un conjunto de vehículos no tripulados

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 Aeroespacial Simulación Dinámica de un Conjunto de Vehículos no Tripulados Autor: Diego Zarco Vaz Tutor: D. Eduardo Fernández Camacho Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Simulación Dinámica de un Conjunto de Vehículos no Tripulados Autor: Diego Zarco Vaz Tutor: D. Eduardo Fernández Camacho Catedrático de Universidad Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Grado: Simulación Dinámica de un Conjunto de Vehículos no Tripulados Autor: Diego Zarco Vaz Tutor: D. Eduardo Fernández Camacho 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 En primer lugar quiero agradecerle a mi tutor D. Eduardo el haber aceptado mi solicitud como alumno de TFG y también que me permitiese tener flexibilidad a la hora de desarrollar el mismo. Después quería agradecerle a mis padres el apoyo que me han dado, y en particular a mi madre que me ha seguido de cerca el desarrollo de este proyecto. Finalmente, estoy muy agradecido con las personas que han estado cerca de mí durante mi camino por este grado, y que me han animado y motivado a seguir avanzando en los momentos más difíciles. Diego Zarco Vaz Sevilla, 2022 I Resumen En este proyecto se desarrollará un algoritmo, utilizando el lenguaje de MATLAB, que diseñe y simule unas trayectorias para un conjunto de vehículos no tripulados. Esta flota de vehículos tendrá dos misiones de inspección que, tras la simulación, serán puestas a prueba en un entorno real. III XÍndice Resultados 63 2.4 Algoritmo completo 64 2.4.1 Función unificadora 64 2.4.2 Entradas a través de un script 69 2.4.3 Entradas a través de MATLAB App Designer 70 3 Simulación 75 3.1 Inspección de un punto 75 3.1.1 Ejemplo 1 77 3.1.2 Ejemplo 2 80 3.1.3 Ejemplo 3 82 3.1.4 Ejemplo 4 86 3.2 Inspección de un área 88 3.2.1 Ejemplo 1 88 3.2.2 Ejemplo 2 92 3.2.3 Ejemplo 3 94 3.2.4 Ejemplo 4 97 4 Pruebas Reales 101 4.1 Introducción 101 4.1.1 Banco de pruebas 102 4.1.2 Crazyflies 102 4.2 1ª Aplicación 103 4.2.1 Simulación 103 4.2.2 Entorno real 104 4.3 2ª Aplicación 105 4.3.1 Simulación 105 4.3.2 Entorno real 106 4.4 3ª Aplicación 107 4.4.1 Simulación 107 4.4.2 Entorno real 108 5 Conclusiones 111 Índice de Figuras 113 Índice de Códigos 117 Bibliografía 119 1 Introducción En este documento, se va a mostrar un nuevo algoritmo que se ha desarrollado para diseñar trayectorias y simularlas. Estas trayectorias tienen como misión recorrer territorios para ejecutar distintas aplicaciones. Estas pueden ser la monitorización de distintas extensiones de terreno, llevar las comunicaciones a sitios remotos o realizar un análisis sobre las zonas inspeccionadas. Para conseguir ejecutar estas misiones, se tienen que imponer una serie de restricciones que las trayectorias tienen que respetar. Además, una vez creadas estas rutas, se va a presentar una simulación en la que una flota de vehículos no tripulados recorrerá los caminos diseñados para comprobar que cumplen con la misión correspondiente. Estas trayectorias y la posterior simulación, se basarán en el uso de vehículos no tripulados conocidos como HAPS (High Altitude Pseudo Satellites). Más adelante se explicará en que consisten estos vehículos. Además, se revisarán algunos planificadores de trayectorias existentes. De esta manera se encontrará el método más conveniente y eficiente para las misiones para las que están destinadas las rutas del presente trabajo. 1.1 Estado del arte En este apartado se van a explicar algunos conceptos necesarios para la realización de este documento. Estos conceptos son los HAPS, ya mencionados, y los planificadores de trayectorias existentes hasta el momento que se consideran más adecuados para el trabajo. 1.1.1 High Altitude Pseudo Satellites En la actualidad existe una zona del espacio por encima de los 20 kilómetros de altitud que no se encuentra aprovechada. En esta zona, también conocida como la estratosfera, es donde aparece el concepto de HAPS (High Altitude Pseudo Satellites) o HAPs (High Altitude Platforms), Figura 1.1. Estos pseudo-satélites son objetos voladores no tripulados que vuelan a grandes altitudes, que tienen como objetivo mantenerse en vuelo durante grandes periodos de tiempo, llegando a estar meses o incluso años en funcionamiento. De esta manera pueden trabajar mucho más tiempo realizando misiones de vigilancia o llevando las comunicaciones a otras zonas del planeta [25]. 1 2Capítulo 1. Introducción Figura 1.1 High Altitude Pseudo Satellites. El concepto de HAPS se estudia desde hace muchos años. Su historia se remonta al siglo pasado, en los 30s ya se comenzaron proyectos que enviaban globos a la estratosfera. Posteriormente, sobre 1950-1960, se experimentó con aeronaves capaces de alcanzar altitudes estratosféricas. Finalmente, a partir de 1990 ya se lanzaban proyectos los cuales tenían como objetivo analizar las posibles aplicaciones de los HAPS en telecomunicaciones y accesos remotos [9]. Detrás de estos proyectos se encuentran muchas empresas famosas las cuales han diseñados sus propios HAPS. Entre ellos están Project Loon de Google, Zephyr de Airbus o Aquila de Facebook. Algunos de estos proyectos se muestran en la Figura 1.2. (a) (b) Figura 1.2 Histórico de proyectos con HAPS. Los HAPS pueden clasificarse en dos tipos. El primer tipo se le conoce como LTA (Lighter Than Air) y corresponde a los vehículos HAPS que son más ligeros que el aire. El segundo tipo es el HTA (Heavier Than Air), que se refiere al vehículo HAPS que es más pesado que el aire. En la Figura 1.3 se muestran ejemplos de ambos tipos de HAPS. Por ejemplo, los globos aerostáticos o los dirigibles son ejemplo de LTA. Mientras que como ejemplo de HTA podrían ser vehículos no tripulados de ala fija. Existen ciertas ventajas de unos tipos respecto a otros pero su uso depende fundamentalmente de la aplicación que se le quiera dar al proyecto. Estas aeronaves se encuentran entre el mundo de los satélites y el de los drones, pudiendo aprovechar las ventajas de ambos mundos y sacando el mayor partido a sus capacidades. Al operar a altas altitudes, unos 20 kilómetros de distancia con la tierra, estos vehículos pueden llegar a cubrir grandes distancias con su campo de visión. Además, no se encuentran tan alejados como los satélites habituales, los cuales están a más de 340 kilómetros. Por tanto, las imágenes realizadas desde los HAPS presentan mayor resolución asemejándose a las obtenidas con un dron. También se puede destacar de este híbrido entre satélite y dron que, al no estar tan alejado de la tierra, no provoca tanto retraso en las comunicaciones como ocurre en los satélites convencionales [21]. 1.1 Estado del arte 3 (a) LTA. Stratobus de Thales. (b) HTA. Zephyr de Airbus. Figura 1.3 Tipos de HAPS. Para que estos pseudo satélites se mantengan mucho tiempo en vuelo necesitan aprovechar la energía solar. A través de paneles solares, recargan sus baterías a lo largo del día y utilizan la energía recolectada en la noche. Además de aprovechar la energía obtenida durante el día, estos vehículos realizan maniobras de ascenso y descenso para optimizar el uso de la misma. Por el día, los HAPS aprovechan la energía del sol para volar a las altitudes deseadas. Posteriormente, cuando llega la noche y el sol se oculta, los HAPS se dejan caer de forma controlada y, de esta manera, usan la menor cantidad de energía posible [ 14 ]. El inconveniente que existe en la actualidad es que se necesita que esta tecnología mejore. Por otro lado, teniendo en cuenta las altitudes a las que vuelan estos vehículos, se observa que en ningún momento interrumpen o causan interferencias con el tráfico aéreo existente en las zonas de vuelo. También, la franja del espacio por la que se desplazan estos vehículos tiene la ventaja de que no se ve afectada por corrientes de viento muy fuertes. Estas corrientes, conocidas como Jet Stream, se encuentran entre los 9 y los 16 kilómetros de altitud [ 8 ] y complicarían el diseño y el funcionamiento de los pseudo-satélites. Además de estas fuertes corrientes de viento, también se evitan gran parte de los bancos de nubes que se encuentran en el cielo. Una vez expuestas las características de este tipo de aeronaves, conviene indicar cuales son las posibles aplicaciones para las que servirían. Las misiones a las que pueden ir destinadas este tipo de aeronaves son muy variadas, Figura 1.4. A continuación se muestra una lista de ejemplos dónde los HAPS podrían tener una amplia repercusión. •Defensa •Vigilancia •Ambiental/Urbano •Navegación •Telecomunicaciones Por ejemplo, cuando se habla de aplicar a los HAPS en misiones de defensa, se refiere a la posibilidad de transmitir comunicaciones o realizar vigilancias de áreas de alta seguridad. Otro ejemplo de la aplicación en defensa podría ser utilizar estos vehículos para controlar las fronteras. Un ejemplo de misiones de vigilancia podría ser el control del tráfico marítimo y aéreo aprovechando las altitudes en las que se trabaja. En otros casos, se podrían utilizar estas aeronaves para el mapeo y la cartografía de ciertas extensiones de terreno. También tiene aplicaciones cuando se habla del ambiente, ya que se pueden realizar trabajos de prevención/vigilancia de incendios o medir la contaminación en ciertas zonas. Por último, en el caso de las telecomunicaciones, se estaría hablando de utilizar los HAPS para dar cobertura en emergencias o catástrofes. 4Capítulo 1. Introducción Figura 1.4 Aplicaciones para las trayectorias [10]. Una de las posibles aplicaciones del presente proyecto podría ser la inspección, monitorización y vigilancia de grandes o pequeños terrenos, ya sean cultivos agrarios o plantas solares. Por ello, se van a utilizar los HAPS de ala fija para el diseño de las trayectorias diseñadas, mientras que la simulación podrá adaptarse a otros drones como los quadrotor. 1.1.2 Planificador de trayectorias En este apartado se realiza un recorrido por los algoritmos planificadores de trayectorias que pueden utilizarse para las aplicaciones que se propondrán en este trabajo. En la actualidad, el término de "Path Planning" ha sido muy utilizado en el mundo de la robótica. Estos algoritmos, de los que hay una gran variedad, buscan una sucesión de posiciones para un robot, que lo lleven desde un estado inicial a un estado final. Estas posiciones dependerán de la distribución de los obstáculos en el entorno. Además, también se tendrá que tener en cuenta la geometría del robot y su maniobrabilidad. Los algoritmos planificadores se suelen diferenciar por la manera de buscar las mejores rutas, por su velocidad a la hora de generar las trayectorias o por la información previa que necesitan para construir esta ruta. Por lo general, lo más adecuado es identificar cuales son las características del problema, y a partir de ahí, seleccionar el método para resolverlo. Los algoritmos de planificación de trayectorias más actuales se caracterizan por obtener la mejor ruta posible o una ruta que se aproxime bastante a la anterior mencionada. Generalmente, la mejor ruta se refiere a la ruta óptima, la cual se consigue minimizando las funciones objetivo del problema. Por ejemplo, se pretende minimizar el tiempo total que se tarda en recorrer la ruta en misiones de búsqueda y rescate. Otro ejemplo podría ser cuando se intenta minimizar la energía que utiliza el vehículo para recorrer la ruta como se puede observar para los casos de vuelos con HAPS. Existen muchas formas de clasificar los algoritmos de planificación. En este trabajo se ha enfocado el estudio en algoritmos cuya planificación no sea dinámica, es decir, no existen obstáculos móviles. Y, además, se conoce el mapa sobre el que se va a trabajar incluyendo los puntos inicial y final. Estos serían los condicionantes para la búsqueda de un tipo de algoritmo de planificación efectivo para el propósito del trabajo. Un tipo de algoritmos que concuerdan con dichas especificaciones son los conocidos como "C Space Search Based Path Planning Algorithms" [ 23 ]. Dentro de ellos se puede distinguir entre "Graph Search" y "Sampling-Based". En la Figura 1.5 se muestra una serie de ejemplos pertenecientes a la categoría de algoritmos "Graph Search". 1.1 Estado del arte 5 Figura 1.5 Algoritmos planificadores. Sin embargo, la categoría "Sampling-Based" parece que se ajusta más a las necesidades del presente trabajo. Dentro de ella se encuentra el algoritmo conocido como RRT (Rapidly Random Tree). Este algoritmo planificador se basa en la construcción de un árbol de puntos que comienzan en el punto inicial y se va expandiendo hasta encontrar el punto final, Figura 1.6. Se le llama random o aleatorio por la forma en la que genera los waypoints por los que pasará el vehículo. Estas posiciones se conocen como nodos y se crean partiendo del nodo inicial. A partir de este, se elige un punto de forma aleatoria en el mapa, se extiende una rama desde el punto inicial en dirección al punto aleatorio y se crea un nodo a una distancia del nodo inicial elegida por el usuario. A continuación, se guarda el nuevo nodo y se sigue expandiendo el árbol de la misma manera hasta alcanzar el punto final. Una vez encontrado el estado final, se realiza una inspección para ver si existe un mejor camino según los nodos existentes. Figura 1.6 Algoritmo RRT. El presente trabajo se centrará en crear un algoritmo que diseñe las trayectorias que van a recorrer las aeronaves desde un punto inicial a un punto final, sin abordar el control del dron durante el vuelo. Los condicionantes de estas trayectorias son, en primer lugar que no existan obstáculos en las mismas, y en segundo lugar, que la flota esté compuesta por aeronaves de ala fija. Esto último implica que las trayectorias no pueden tener giros bruscos ya que en este tipo de vehículos se necesita cierto espacio para poder maniobrar manteniendo la seguridad. Debido a estas dos premisas, los algoritmos planificadores mencionados en este apartado, aun aportando mucha información, no se ajustan a las necesidades del proyecto. En un principio se comenzó creando trayectorias utilizando una variante del algoritmo RRT combinado con un tipo de curvas denominadas Caminos de Dubins. Pero dado la ausencia de obstáculos fijos o móviles, se observó que no era necesario recurrir a algoritmo de planificación que busquen el camino óptimo basándose en la distribución de los obstáculos. En el caso que se plantea en este trabajo, se le da más importancia a la maniobrabilidad de la flota de vehículos no tripulados y al tipo de elementos que se quieren visualizar. Por todo ello se decidió crear un algoritmo diferente, de creación propia. 6Capítulo 1. Introducción 1.2 Objetivos En este proyecto se desarrollará un algoritmo, utilizando el lenguaje de MATLAB, que diseñe y simule unas trayectorias para un conjunto de vehículos no tripulados. Esta flota de vehículos tendrá dos misiones de inspección que deberán llevar a cabo de forma eficiente y evitando colisiones entre las distintas aeronaves. La primera misión se corresponde con la vigilancia de un punto del mapa considerando algunas restricciones como las características de la aeronave o el número de inspecciones que se quieran realizar. La segunda misión consiste en la visualización de una zona grande de terreno durante un tiempo determinado. Finalmente, tras obtener las trayectorias se realiza una simulación en la que un conjunto de vehículos no tripulados recorren estas rutas para verificar que se cumplen todas las condiciones impuestas por el usuario. En relación a los objetivos a inspeccionar, puede tratarse desde una planta solar de 100 hectáreas, como podría ser Solnova en Sanlúcar la Mayor, sobre la que se quiere realizar una inspección, hasta una comunidad autónoma entera como es Andalucía. En la Figura 1.7 se muestran los ejemplos mencionados hacia los que se podrían orientar las misiones. (a) Planta solar Solnova. (b) Mapa de Andalucía. Figura 1.7 Posibles áreas a inspeccionar. Debido a la variabilidad en las dimensiones del terreno a inspeccionar, el algoritmo tiene que ser escalable, es decir, tiene que crear trayectorias de manera adecuada independientemente de si la misión le concierne a un HAPS que vuela a muchos kilómetros de altitud o si la misión la ejecuta una flota de vehículos no tripulados sobre una planta solar. Finalmente, la simulación correspondiente tiene que aproximarse lo más posible a la realidad. Será necesario comprobar que no existen colisiones entre las aeronaves, ya que aunque no haya obstáculos en las condiciones del problema, si hay que tener en cuenta que se va a utilizar una flota de vehículos los cuales pueden colisionar entre sí. Por tanto, las trayectorias creadas tendrán componentes que permitan garantizar la seguridad de los vehículos. Al final del trabajo, el funcionamiento del algoritmo se comprobará realizando pruebas reales con las trayectorias creadas. 1.3 Estructura del documento 7 1.3 Estructura del documento Una vez introducidos, en el capítulo uno, los conceptos necesarios para el desarrollo del presente trabajo y establecidos los objetivos, a continuación se indica la estructura que tendrá el resto del documento, la cual corresponde a tres capítulos y unas conclusiones. El segundo capítulo de este documento consistirá en el desarrollo del algoritmo, para lo que será necesario explicar dos conceptos: las curvas de Dubins, que se utilizarán como base para la creación de las trayectorias, y el campo de visión de una aeronave, importante para determinar la forma de las rutas diseñadas. Posteriormente, se hará una breve introducción al primer caso de uso donde se concretará la misión correspondiente y se expondrán las dificultades que se fueron planteando y que moldearon la forma de su desarrollo. A continuación, se desarrollará el algoritmo correspondiente a la inspección de un punto. El mismo esquema se seguirá para la exposición del segundo caso de uso, la inspección de un área. Por último, se creará un único código que unificará ambos algoritmos y una aplicación donde, de forma intuitiva, se introducen los datos de partida para la ejecución de la misión. En el tercer capítulo se explicará el código correspondiente a la simulación de las trayectorias y se expondrán cuatro ejemplos para cada caso de uso de distintas situaciones que muestran el potencial del algoritmo. El cuarto capítulo se corresponde con las pruebas reales. Se comienza describiendo el banco de pruebas y los drones crazyflies utilizados en los ensayos. Posteriormente, se presentan tres aplicaciones. Para ello, primero se define el objetivo, después se crea la trayectoria y se simula, y por último, se verifican los resultados de la simulación con una prueba en un entorno real. Finalmente, se incluirán las conclusiones y posibles líneas futuras de trabajo. 2 Algoritmo En esta sección primero se hablará de dos conceptos muy importantes de cara al problema que se pretende resolver. Estos conceptos son las curvas de Dubins y el campo de visión de una aeronave. Además de explicar a fondo estas dos ideas, se comentará como se pretenden incluir en el trabajo que se desarrollará en el resto del capítulo. Posteriormente, se procederá a introducir primero el caso de la inspección de un punto. En este problema, se tiene como objetivo vigilar un cierto punto en el mapa bajo unas condiciones de velocidad y tiempo. Por otro lado, también se explicará el segundo problema planteado en este trabajo, la inspección de un área. Al igual que en el primer caso, bajo unas restricciones de velocidad y de tiempo, se vigilará un área extensa sobre el mapa con un conjunto de vehículos no tripulados. 2.1 Conceptos previos 2.1.1 Curvas de Dubins El concepto de curva de Dubins o camino de Dubins (en inglés, Dubins curve o Dubins path) se refiere a la curva de longitud mínima con una restricción en la curvatura y con tangentes y posiciones inicial y final prescritas, considerando que el vehículo que recorre la ruta solo puede ir hacia delante. Si el vehículo pudiera ir hacia atrás, estaríamos hablando de las curvas Reeds-Shepp, tema que no se tratará en este documento [ 5 ]. En la Figura 2.1 se puede observar un ejemplo simple de una curva de Dubins. Más adelante, se explicarán los elementos que la componen. Figura 2.1 Curva de Dubins. Ejemplo básico. 9 16 Capítulo 2. Algoritmo Figura 2.12 Cálculo curva de Dubins RSL (6). Finalmente, el camino de Dubins se compone de tres tramos. El primer tramo se encuentra en la circunferencia R1 , sería el arco que se describe desde el punto inicial  S hasta el punto de tangencia  T1 , recorrido en sentido horario. El segundo tramo sería el segmento  RSL ya calculado, que corresponde con la tangente interna a ambas circunferencias que une los puntos  T1 y  T2 . Finalmente, el último tramo es el trozo de circunferencia que se recorre, en sentido antihorario, desde el punto de tangencia  T2, hasta el punto final  E. En la Figura 2.13 se muestra como quedaría la trayectoria final. Figura 2.13 Cálculo curva de Dubins RSL (7). De forma análoga, se podrían obtener las cinco curvas de Dubins restantes para los mismos puntos inicial y final. Los cálculos varían respecto a este ejemplo, pero la metodología es similar. Cabe destacar que no siempre será posible calcular los seis caminos de Dubins. Ya que en los casos formados por tres circunferencias consecutivas, RLR y LRL, solo es posible el cálculo de estas curvas si la distancia entre los puntos inicial y final es menor de seis veces el radio mínimo de giro. Esto se debe a que es la máxima distancia alcanzable por tres circunferencias tangentes entre sí. Función de MATLAB Una vez explicado el modo de calculo de los caminos de Dubins, se prosigue mostrando cómo se podrían calcular estas curvas en el software de MATLAB. Todos los algoritmos desarrollados a lo largo de este proyecto, se han enfocado en el lenguaje de MATLAB, aunque estos tienen adaptabilidad para otros lenguajes como se comentará más adelante. En MATLAB, se puede utilizar una función específica que te encuentra las curvas de Dubins dados un punto inicial y un punto final, pudiendo cambiar parámetros como el radio mínimo de giro. El código necesario para ejecutar esto se puede encontrar en la propia página de MATLAB [ 6 ] y corresponde al Código 2.1. Código 2.1 Función de curva de Dubins de MATLAB. 1% Create a dubinsConnection object. 2.1 Conceptos previos 17 2dubConnObj = dubinsConnection; 3 4% Define start and goal poses as [x y theta] vectors. 5startPose = [0 0 0]; 6goalPose = [1 1 pi]; 7 8% Calculate a valid path segment to connect the poses. 9[pathSegObj, pathCosts] = connect(dubConnObj,startPose,goalPose); 10 11 % Show the generated path. 12 show(pathSegObj{1}) En este ejemplo se pretende obtener el camino de Dubins más corto entre el punto inicial, el cuál corresponde al [0,0,0] , es decir, coordenadas x=0 , y=0 y con un ángulo de salida de 0◦ respecto a la horizontal; y un punto final [1,1,π] cuyo ángulo final tiene que ser de 180◦ respecto a la horizontal. Para obtener la curva de Dubins primero hay que crear un objeto que MATLAB denomina dubinsConnection como se muestra en la línea 2 del Código 2.1). Posteriormente, se definen los puntos inicial y final con el formato ya indicado [x,y,θ] . Por último, se hace uso de la función connect, línea 9 del Código 2.1, que conecta los puntos siguiendo un camino de Dubins. Cabe destacar que, como parámetros predeterminados se encuentra por ejemplo el radio mínimo de giro con una unidad de MATLAB como valor por defecto. Finalmente, en la Figura 2.14 se obtiene una visualización de los resultados obtenidos gracias al uso del comando show de MATLAB. Figura 2.14 Ejemplo en MATLAB de curva de Dubins. De la gráfica obtenida se pueden identificar varios elementos. El punto inicial corresponde al objeto situado en el [0,0] , pintado de color verde. El punto final por otro lado, se encuentra en la posición [1,1] de la figura, pintado en rojo. Otro aspecto de la leyenda a destacar es la ruta o camino que se encuentra como un trazo de color naranja etiquetado en la leyenda como path. A lo largo de este trazo existen dos puntos negros con una pequeña línea negra. Estos elementos indican la posición de los puntos de tangencia y la dirección que seguiría el vehículo al pasar por ese punto. Se puede observar también que la curva calculada sería del tipo RLR, es decir, una primera curva en sentido horario, posteriormente se recorre un trozo de circunferencia en sentido antihorario para acabar con otro pequeño trozo en sentido horario. La función de MATLAB te indica que tipo de curva es, la longitud de la curva, el coste del 18 Capítulo 2. Algoritmo camino y más información que se guarda en las variables de salida de la línea 9 del Código 2.1. El elemento pathSegObj consiste un una celda que incluye los puntos inicial y final, el radio mínimo de giro, la longitud de cada uno de los tres tramos, el tipo de curva de Dubins y por último la longitud total de la trayectoria. A continuación, se van a realizar unas pequeñas modificaciones en el código para poder obtener un ejemplo dónde se obtengan todas las curvas posibles. Código 2.2 Función de curva de Dubins de MATLAB modificada. 1% Create a dubinsConnection object 2dubConnObj = dubinsConnection; 3 4% Change minimum turning radius 5dubConnObj.MinTurningRadius = 0.5; 6 7% Define start and goal poses as [x y theta] vectors 8startPose = [0 0 0]; 9goalPose = [1 0 -pi/4]; 10 11 % Calculate a valid path segment to connect the poses 12 [pathSegObj, pathCosts] = connect(dubConnObj,startPose,goalPose," PathSegments","all"); 13 14 % Show the generated path 15 for i = 1:6 16 figure() 17 show(pathSegObj{i}) 18 axis equal 19 end Los aspectos más importantes que se han modificado para obtener el Código 2.2 son los siguientes. Se ha cambiado el punto final a la posición [1,0,−π 4] y se ha disminuido el radio mínimo de giro a la mitad, por tanto, tendrá un valor de 0.5 unidades. Por último, se han modificado las entradas de la función connect para obtener todas las curvas y se ha añadido un bucle que permita sacar una figura para cada camino de Dubins. La Figura 2.15 y la Figura 2.16 muestran los resultados obtenidos al ejecutar el Código 2.2. (a) RSR. (b) RSL. (c) LSR. Figura 2.15 Ejemplo modificado en MATLAB de curvas de Dubins (1). 2.1 Conceptos previos 19 (a) LSL. (b) RLR. (c) LRL. Figura 2.16 Ejemplo modificado en MATLAB de curvas de Dubins (2). Todas las imágenes comparten el mismo punto inicial, punto final y radio mínimo de giro. Los elementos que aparecen en las figuras y en la leyenda de cada uno son los mismos elementos que se obtuvieron en la Figura 2.14. Se puede observar que hay caminos de Dubins que no parecen óptimos comparados con otros. Póngase un ejemplo: en la Figura 2.15a la ruta realiza casi una vuelta completa en sentido horario para alcanzar el ángulo necesario que le permite aproximarse al punto final. Por otro lado, en la Figura 2.15c se observa como sería posible llegar desde el mismo punto inicial hasta el final, utilizando 3 tramos muy cortos y directos hacia la meta deseada. Es por ello, que una curva de Dubins cualquiera no siempre será el camino más corto, sino que hay que tener en cuenta cuál de las seis posibilidades es la más corta. Otros ejemplos a destacar en este caso, serían los caminos de Dubins compuestos solo por curvas, y por tanto, no constan de ningún tramo rectilíneo. En los casos de la Figura 2.16b y de la Figura 2.16c se observa que se completa la ruta con tres tramos curvilíneos formados por arcos de circunferencias de radio mínimo. Para el ejemplo propuesto, esta manera de actuar no es la adecuado ya que te obliga a recorrer gran parte de las circunferencias para hallar el ángulo adecuado para pasar al siguiente tramo. A pesar de la efectividad de la función de MATLAB, es difícil tener acceso al cálculo de estas trayectorias de una forma más profunda. Es por ello, que se planteó la idea de buscar unas funciones sobre las cuales se tenga más control del propio algoritmo. De esta manera, se podrán hacer cambios a la función según los intereses de la investigación. Por tanto, se buscó una función que al igual que la de MATLAB, calculara los caminos de Dubins con la restricción en el radio de giro de la curva. Función de Andrew Walker En 2008, Andrew Walker, realizó un software que calcula el camino más corto para un vehículo de Dubins, es decir, que tiene una restricción de maniobrabilidad y solo va hacia delante. Andrew comenta que sus investigaciones se basan en los trabajos de [ 15 ]y[ 20 ], y que el enfoque usado para el cálculo de estas rutas sería el encontrado en el artículo [24] del cual se hablará más adelante. El código que realizó Andrew Walker fue utilizando el lenguaje C y lo publicó en la web de GitHub [ 26 ]. Esta web es un lugar en internet creado para alojar diferentes proyectos a los cuales se puede acceder y seguir investigando con ellos. Más adelante, en 2016 Hsin-Yi Kang o también conocido en la plataforma como EwingKang, decidió traducir el código de Andrew al lenguaje de MATLAB [ 17 ]. Además, Hsin-Yi tiene varios proyectos subidos a la web hablando sobre las curvas de Dubins y algoritmos de planificación de trayectorias como el RRT del que se habló con anterioridad. Finalmente, en este trabajo se va a analizar la traducción que hizo Ewingkang en MATLAB del código que creó Andrew Walker sobre las curvas de Dubins. Entre los archivos que se encuentran en la publicación de Hsin-Yi encontramos dos documentos con formato ".m" que servirán como núcleo del desarrollo de las trayectorias en este proyecto. Estos 20 Capítulo 2. Algoritmo archivos serían dubins_curve.m y dubins_core.m . Primero, se van a introducir las funciones en un script de MATLAB tal como recomienda el autor. Posteriormente se ejecutarán para ver los resultados obtenidos. En la Figura 2.17 se muestra los resultados correspondientes al ejecutar el Código 2.3. Código 2.3 Ejemplo del código de Hsin-Yi. 1P1 = [1 1 0]; 2P2 = [4 4 -pi/4]; 3r = 1; 4stepsize = 0.01; 5PATH = dubins_curve(P1, P2, r, stepsize); 6axis([0 5 0 5]), daspect([1 1 1]) Figura 2.17 Ejemplo del código de Hsin-Yi. Para utilizar las funciones primero hay que crear los puntos inicial y final de la trayectoria. En la Figura 2.17 están marcados con las palabras start para el punto inicial y end para el punto final. Después, se le indica a la función cual es el radio mínimo de giro que se use como restricción para construir el camino. Y, por último, hay que indicar como entrada de la función el "stepsize" o la distancia entre puntos que se quiere a la hora de realizar la trayectoria. Esta distancia sirve para obtener una mayor precisión en la curva a causa de utilizar mayor capacidad de cómputo. Además de lo indicado, se ha añadido una línea final (línea 6 del Código 2.3) dónde se indican las dimensiones de la figura de salida y la relación de aspecto entre los distintos ejes. Esto sirve para que las curvas se vean como trayectorias circulares. Finalmente, se obtiene como salida de la función una matriz de 3 columnas y, en este caso, 525 filas que corresponderían a los distintos puntos de la trayectoria. Las tres columnas son la coordenada en x , la coordenada en y y el ángulo θ , en radianes, respectivamente. Además, se muestra en la ventana de comandos de MATLAB el mensaje que se ve en la Figura 2.18. Figura 2.18 Mensaje del código de Hsin-Yi. 2.1 Conceptos previos 21 En este mensaje se muestra el tiempo que se ha tardado en computar la trayectoria y el tiempo que se ha tardado en mostrar la figura correspondiente. La función ejecutada, tiene un parámetro adicional de entrada llamado quiet. Este parámetro se encuentra de forma predeterminada con un valor false o falso dentro de la función. Si se le introduce un valor de true o verdadero en las entradas, esto produce que no se muestre por pantalla un dibujo de como sería la trayectoria. Esta idea será utilizada en posteriores desarrollos. A continuación, se procede a analizar la función que se ha ejecutado: dubins_curve.m . Ya que, para poder utilizarla en futuras investigaciones, lo más conveniente es entender su funcionamiento a la perfección. Al abrir la función en MATLAB, las primeras líneas corresponden a una explicación de como utilizar la función. Además, se añaden las referencias a la web GitHub y más datos sobre el autor. Código 2.4 Función dubins_curve.m (1). 1%DUBINS_CURVE Find the Dubins path (shortest curve) between two points. 2% PATH = DUBINS_CURVE(P1, P2, r, stepsize) finds the shortest curve 3% that connects two points in the Euclidean plane with a constraint of 4% the curvature of the path. The start and finish orientations P1 and 5% P2 are defined as [x, y, theta]. The turning radius (r) and stepsize 6% will be defined automatically if their value is <= 0. The output 7% PATH is an [mx3] array consisting of m rows of [x, y, theta] values. 8% 9% PATH = DUBINS_CURVE(P1, P2, r, stepsize, quiet) performs the same as 10 % above, however if quiet == true, then no plots of command window 11 % output will be generated. Ommitting this input will result in quiet 12 % = false/0. 13 % 14 % This function handles the interface to dubins_core.m to give a more 15 % intuitive tool for finding the Dubins path. Como se observa en el Código 2.4, el autor explica cómo usar la función y los parámetros de entrada que la componen. Después, indica como se usa el último parámetro de entrada a la función que, como ya se ha dicho, este sirve para sacar un gráfico con el dibujo de la trayectoria. Finalmente, comenta que dentro de esa función se utiliza el otro archivo hallado en la librería, dubins_core , como herramienta para encontrar el camino de Dubins. Tras estos comentarios iniciales comienza el código correspondiente a la función. De este, solo se hablará de las partes más relevantes del mismo, pasando por encima algunas líneas menos importantes. La primera parte importante se encuentra en las líneas presentadas en el Código 2.5. Código 2.5 Función dubins_curve.m (2). 1% main function 2param = dubins_core(p1, p2, r); 3if stepsize <= 0 4stepsize = dubins_length(param)/1000; 5end 6path = dubins_path_sample_many(param, stepsize); En estas líneas se hace una llamada a la función dubins_core en la cuál se realizan los cálculos 22 Capítulo 2. Algoritmo de los diferentes tramos del camino de Dubins más corto. Estos tramos, salen en forma de valores numéricos en la estructura denominada param , la cuál se explicará más adelante. Tras obtener los parámetros, se analiza el valor dado por el usuario de stepsize y se ejecuta la función interna en el algoritmo denominada dubins_path_sample_many . Esta función convierte los parámetros calculados en el código de dubins_core en puntos a lo largo de la trayectoria separados una distancia igual al stepsize dado. El resto de la función de dubins_curve consiste en realizar el cálculo de las trayectorias ya comentadas a partir de los parámetros y dibujar la gráfica correspondiente. Es por ello, que a continuación se va a comentar como se han llegado a obtener los tramos de los caminos de Dubins dentro de la función dubins_core.m. Al abrir este nuevo archivo, aparecen los mismos comentarios encontrados en la función anterior. Posteriormente, se explican como se han definido los diferentes tramos posibles para las curvas de Dubins. Código 2.6 Función dubins_core.m (1). 1%%%%%%%%%%%%%%%%%% DEFINE %%%%%%%%%%%%%%%%% 2% Here are some usefuldefine headers for better implementation 3% there are 6 types of dubin’s curve, only one will have minimum cost 4% LSL = 1; 5% LSR = 2; 6% RSL = 3; 7% RSR = 4; 8% RLR = 5; 9% LRL = 6; 10 11 % The three segment types a path can be made up of 12 % L_SEG = 1; 13 % S_SEG = 2; 14 % R_SEG = 3; 15 16 % The segment types for each of the Path types 17 %{ 18 DIRDATA = [ L_SEG, S_SEG, L_SEG ;... 19 L_SEG, S_SEG, R_SEG ;... 20 R_SEG, S_SEG, L_SEG ;... 21 R_SEG, S_SEG, R_SEG ;... 22 R_SEG, L_SEG, R_SEG ;... 23 L_SEG, R_SEG, L_SEG ]; 24 %} 25 %%%%%%%%%%%%%%%% END DEFINE %%%%%%%%%%%%%%%% Como se comenta en el Código 2.6, existen seis posibles caminos de Dubins, de los cuales solo uno será el más corto. El sistema para obtener el camino más corto se basará en asignarle un coste a los tramos de cada camino según la longitud de los mismos. Al igual que con el archivo de dubins_curve.m , en este caso se van a obviar las partes menos importantes. A continuación, en el Código 2.7 se muestran el cálculo de todas las posibles curvas de Dubins y entre ellas se busca la más corta. 2.1 Conceptos previos 23 Código 2.7 Función dubins_core.m (2). 1% Second, we find all possible curves 2best_word = -1; 3best_cost = -1; 4test_param(1,:) = dubins_LSL(alpha, beta, d); 5test_param(2,:) = dubins_LSR(alpha, beta, d); 6test_param(3,:) = dubins_RSL(alpha, beta, d); 7test_param(4,:) = dubins_RSR(alpha, beta, d); 8test_param(5,:) = dubins_RLR(alpha, beta, d); 9test_param(6,:) = dubins_LRL(alpha, beta, d); 10 11 for i = 1:1:6 12 if(test_param(i,1) ~= -1) 13 cost = sum(test_param(i,:)); 14 if(cost < best_cost) || (best_cost == -1) 15 best_word = i; 16 best_cost = cost; 17 param.seg_param = test_param(i,:); 18 param.type = i; 19 end 20 end 21 end En este Código 2.7 se calculan individualmente todos los posibles caminos de Dubins. Los parámetros de entrada corresponden a ángulos calculados a partir de los puntos inicial y final. Estos ángulos denominados alpha y beta aparecen explicados en el artículo [ 24 ] que sirvió como guía para la creación de este archivo dubins_core . En el artículo mencionado, se muestra el ejemplo que se va a incluir a continuación. En el ejemplo que se muestra en el Código 2.8, se va a calcular la longitud de los diferentes tramos del camino de Dubins LSL. Código 2.8 Función dubins_core.m (3). 1function param = dubins_LSL(alpha, beta, d) 2tmp0 = d + sin(alpha) - sin(beta); 3p_squared = 2 + (d*d) -(2*cos(alpha - beta)) + (2*d*(sin(alpha) - sin(beta))); 4if( p_squared < 0 ) 5param = [-1, -1, -1]; 6return; 7else 8tmp1 = atan2( (cos(beta)-cos(alpha)), tmp0 ); 9t = mod((-alpha + tmp1 ), 2*pi); 10 p = sqrt( p_squared ); 11 q = mod((beta - tmp1 ), 2*pi); 12 param(1) = t; 13 param(2) = p; 14 param(3) = q; 15 return ; 16 end 24 Capítulo 2. Algoritmo 17 end En este código, a partir de los parámetros alpha , beta y d , se calcula la longitud de los diferentes tramos que componen el camino de Dubins LSL. Estos tramos corresponden a un tramo curvilíneo recorrido en sentido antihorario, un tramo recto y otro tramo igual que el primero. Estas longitudes corresponden en el código a los parámetros t , p y q . Al igual que se hace con este ejemplo, unos cálculos parecidos se realizan para el resto de curvas de Dubins. Todos estos cálculos se encuentran indicados en [24], el cuál ha sido la base para el desarrollo de estos algoritmos. Finalmente, se da por concluida la explicación de los algoritmos que servirán como núcleo en el desarrollo futuro de este trabajo. A continuación se explicará otro concepto muy importante que se utilizará a lo largo del proyecto, el campo de visión de una aeronave. 2.1.2 Campo de visión El campo de visión de una persona es el rango que una persona puede ver a través de los ojos. Es decir, es donde se aprecian los objetos situados en el espacio. En las personas, el rango de visión ronda los 180◦ aunque el rango de visión efectivo o el lugar dónde las personas distinguen bien los elementos no es superior a los 30◦ [ 3 ]. En la Figura 2.19 se muestra el campo de visión de una persona. Figura 2.19 Campo de visión de una persona. Al campo de visión de forma general también se le conoce por su traducción al inglés, Field of View (FOV). Esta visión consta de dos componentes, el campo de visión horizontal y el campo de visión vertical formando ambas un rectángulo en su conjunto [ 27 ]. Al igual que las personas, existen sistemas como las cámaras que tienen su propio campo de visión. Estos sistemas utilizan distintos objetivos según el uso que se les quiera dar. Por ejemplo, se necesitará un objetivo tipo gran angular para campos de visión grandes. Además, esto también suele perjudicar lo que se conoce como la resolución. Es lógico pensar que cuanto más grande es el campo de visión de un objeto, menor es la resolución del mismo, es decir, las pequeñas cosas que aparecen en la imagen serán menos precisas [ 13 ]. Esta precisión de la cámara o de los sensores utilizados aporta las medidas más pequeñas que pueden representarse en un píxel. A lo largo del trabajo se tratará con parámetros mucho más generales que los mencionados aquí. De estos parámetros se hablará a continuación. En misiones de vigilancia medioambiental, de búsqueda y rescate, o muchas otras de las que se habló cuando se mencionó el concepto de HAPS (High Altitude Pseudo Satellites), era sencillo 2.1 Conceptos previos 25 obtener un amplio campo de visión debido a la gran altitud a la que vuelan estas aeronaves. Además, es interesante tener controlado el FOV para saber en todo momento las distancias que abarcan los sensores de estos vehículos no tripulados. A continuación se va a mostrar un esquema en la Figura 2.20 sobre el que se basará el campo de visión utilizado para realizar las trayectorias a lo largo de este proyecto. Figura 2.20 Campo de visión de una aeronave. Esta imagen fue sacada de [ 19 ]. Este artículo titulado "Automated UAV tasks for search and surveillance" será muy recurrente a lo largo del trabajo ya que menciona términos y propone desarrollos muy interesantes a tener en cuenta. En la Figura 2.20 se observa cual es el campo de visión de una aeronave que vuela a una altitud halt con dos perspectivas distintas. En la imagen de la derecha podemos observar una vista de perfil, donde se muestra el rango de visión que tiene la nave hacia delante y hacia detrás. En la otra vista, imagen del lado izquierdo, se observa una perspectiva desde arriba donde se ve el ancho del campo de visión total. Para obtener las dimensiones totales del campo de visión, a parte de la altitud, se necesitan tres ángulos representados en estos esquemas. Estos ángulos son conocidos como el ángulo de elevación del gimbal del sensor ( θg ), el ángulo vertical del FOV ( ηv ) y, por último, el ángulo horizontal del FOV ( ηh ). A partir de estos ángulos y de la altitud de la aeronave, se pueden obtener los parámetros que interesan para calcular las dimensiones sobre el mapa del campo de visión como se detalla en las siguientes ecuaciones. dl=halt tan(θg−ηv/2)(2.11) dt=halt tan(θg+ηv/2)(2.12) dw=2halt cos(θg)tan(ηh/2)(2.13) Como se puede observar, el cálculo de estos parámetros es puramente geométrico. Estos parámetros de interés son la longitud hasta el borde frontal del campo de visión ( dl ), también conocido como leading edge; la longitud hasta el extremo posterior del campo de visión ( dt ) o también trailing edge y, por último, el ancho del campo de visión o width (dw). En los próximos desarrollos, se considerará que el campo de visión tendrá geometría de rectángulo 32 Capítulo 2. Algoritmo caso de utilizar una forma de infinito en la trayectoria, se dedujo que se necesitaría otra aeronave en el caso en el que el radio de giro correspondiente al tiempo impuesto, sea menor que el radio de giro mínimo. Cuando ocurra esto, será necesario introducir otra aeronave en la trayectoria. Posteriormente, según se desarrollaba el código, se vio conveniente descartar esta idea del patrón en forma de infinito, y volver a la forma de hipódromo. El motivo de esto es que era mucho más sencillo controlar los tiempos con las distancias que se alargaban o acortaban en el patrón de espera. Esto ocurría, ya que eran tramos rectilíneos y no consistían en curvas dónde calcular la longitud que se recorría no era tan trivial. Además, el patrón de infinito no estaba compuesto por tramos que formasen circunferencias perfectas, sino que combinaba tramos de circunferencia con tramos rectilíneos. Asimismo, también se podía cumplir el último de los objetivos. El objetivo de identificar cuántas aeronaves eran necesarias para cumplir las especificaciones de tiempo. Esto era posible ya que, en vez de medirlo con el radio de giro mínimo, se medía con la longitud del tramo de aproximación. Se recuerda, que el tiempo por vuelta se controlaba alargando o acortando el tramo de aproximación al punto de inspección. Si este tramo, no se pudiera acortar más porque ese valor se hiciera negativo, significaría que se necesita un nuevo UAV para cumplir la condición de tiempo. En la Figura 2.27, se puede observar como quedaría una trayectoria donde el tramo de aproximación sea mínimo y si se disminuyera el tiempo de inspección, habría que añadir otra aeronave. Figura 2.27 Inspección con espera en hipódromo y tiempo mínimo. Finalmente, estos fueron todos los objetivos que se fueron proponiendo y que le daban forma a la resolución del problema planteado. A continuación, se va a explicar el código completo donde quedan reflejado todas las cosas de las que se ha estado hablando. 2.2.2 Desarrollo en MATLAB Argumentos de entrada En esta sección del documento se va a comentar el código completo que calcula las trayectorias correspondiente a la inspección de un punto. Primero se van a comentar cuales son las variables que se introducen en la función como argumentos de entrada. Además, se van a indicar los valores que pueden y no pueden tomar. Estos datos son los siguientes: • Estratopuertos. Como ya se ha explicado, los estratopuertos son instalaciones con pistas adecuadas para el despegue y aterrizaje de HAPS (High Altitude Pseudo-Satellite). Este dato se introducirá en el algoritmo como una serie de puntos de 3 coordenadas. 2.2 Caso de uso 1: Inspección de un punto 33 Las dos primeras coordenadas corresponderán con su posición en el plano xy . Estas coordenadas pueden ser indicadas en cualquier unidad de distancia mientras sea coherente con el resto de datos, aunque lo habitual es utilizar metros. La tercera coordenada será la orientación de la aeronave respecto al eje x en forma de ángulo en grados. Se irán incluyendo estos vectores de coordenadas en una matriz. Esta matriz será el dato a proporcionar como entrada al algoritmo. Se pueden incluir todos los estratopuertos que el usuario quiera. Aunque, el algoritmo hará uso solo de los necesarios según el número de UAVs que vayan a volar. Además, si el número de estratopuertos es menor que el de aeronaves necesitadas, se volverán a usar los estratopuertos en el orden en el que se hayan introducido. • Punto de inspección. En este primer caso de uso, se considerará un solo punto con el mismo formato que las coordenadas de los estratopuertos [x,y,θ] . Sobre este punto se realizará la inspección. También cabe destacar que la orientación que se introduzca en el punto de inspección servirá como dirección del circuito de espera. Al igual que en los estratopuertos, las coordenadas x e y se deben introducir con las mismas unidades que las del punto anterior, y la orientación en radianes. • Velocidad de la aeronave. Como su propio nombre indica, el tercero de los argumentos de entrada de la función corresponde con la velocidad que llevará la aeronave en vuelo. Este valor se considerará constante a lo largo de la trayectoria. Las unidades de esta velocidad deben corresponder con las unidades de distancia utilizadas anteriormente y con las unidades de tiempo que quieran utilizarse en el futuro. Lo habitual es que la velocidad de la aeronave sea metros por segundo. Los valores típicos de los HAPS rondan los 30 m/s y nunca se admitirán valores negativos para la velocidad. En caso de que se introduzca un valor negativo, el código mostrará por pantalla un mensaje de error indicando el origen del mismo. • Radio de giro mínimo de la aeronave. Este radio será el dato que se introducirá en las curvas de Dubins para que las trayectorias respeten este valor al realizar giros y curvas. Este valor no tiene que ser el valor mínimo que puede alcanzar la aeronave al realizar un giro, ya sea por fallo estructural o para evitar entradas en perdidas durante el vuelo. Ya que, también se puede configurar para que este valor represente el radio que se quiere seguir a la hora de realizar un viraje. Al igual que con la velocidad de la aeronave, este valor no puede tomar valores negativos ya que no tendría un sentido lógico. Es por ello, que si el usuario introduce un valor que no sea positivo, el algoritmo lanzará un mensaje con un error comentando que esto no es posible. • Tiempo total de vuelo. Este valor de tiempo corresponde con el tiempo que se quiere que la aeronave esté en el aire. Realizando una serie de cálculos que se verán más adelante, se obtiene el número de vueltas al circuito necesarias que representan este tiempo total de vuelo. Para el cálculo de este número de vueltas no se tiene en cuenta el tiempo que la aeronave tarda en llegar o en volver al estratopuerto, únicamente el vuelo en el circuito de inspección. Este valor temporal obviamente no puede tener valores menores o iguales a 0. Por lo tanto, en cuanto el código detecta que se ha introducido un valor que no es el adecuado, al igual que en los puntos anteriores, mostrará un mensaje indicando que hay que cambiar el valor del tiempo total debido a que no es adecuado para el caso de uso. • Tiempo de inspección. Este valor de tiempo, a diferencia del anterior, corresponde con la frecuencia con la que se quiere que la aeronave pase por el punto a inspeccionar. Es decir, el tiempo que se quiere que la aeronave tarde en dar una vuelta completa al circuito. Sabiendo el tiempo de vuelo total y el tiempo de vuelta, se puede obtener el número de vueltas 34 Capítulo 2. Algoritmo que será necesario dar para respetar el tiempo total de vuelo. Estos dos valores temporales se tienen que dar con las unidades correspondientes a las que se utilizaron con la velocidad, normalmente será en segundos. En esta ocasión vuelve a ocurrir que un valor 0 o negativo de tiempo no es adecuado y, por tanto, se mandará un mensaje de error al usuario. • Campo de visión. El siguiente de los argumentos de entrada es el campo de visión o Field of View. Los parámetros que se iban a usar ya se explicaron en el apartado anterior donde se hablaba sobre el campo de visión de una aeronave. Este dato se introducirá como un vector de 3 componentes, donde la primera componente corresponde con el ancho del campo de visión ( dw ). La segunda de las componente sería la longitud del extremo posterior del campo de visión ( dt ). Por último, la tercera componente del vector es la longitud del extremo frontal del campo de visión (dl). Como ya se comentó, la longitud del borde inferior puede tener valores tanto positivos como negativos, a diferencia de el ancho del campo de visión o de la longitud del borde superior que solo pueden tomar valores superiores a 0. En los casos en los que no sea así, al igual que en los argumentos anteriores, se lanzará un mensaje de error con las indicaciones correspondientes. • Altitud. La altitud de la aeronave se indica en este momento y se considera constante durante todo el vuelo. Esto es debido a que estas trayectorias se centran simplemente en la ruta que se realiza sobre el plano. Aun así, se introduce la altitud de la aeronave la cual fue necesaria en pruebas reales que se mostrarán más adelante. Sobre la altitud también se comprueba que tenga valores positivo ya que no tendría sentido dar una coordenada negativa en el eje z que indique la altitud de la aeronave. Por tanto, se comprueba este dato y se emite un mensaje de error en caso de que no se cumpla el requerimiento. • Stepsize. Como ya se comentó en el desarrollo de los algoritmos correspondientes a las curvas de Dubins, el stepsize se denomina como la distancia que existe entre punto y punto de la trayectoria. Cuando menor sea el valor del mismo mayor será la precisión de la trayectoria, esto se ve reflejado sobretodo en las curvas. Aunque un valor del stepsize muy pequeño hace que los resultados tarden más en obtenerse y que el número de puntos en la trayectoria crezca mucho. Este valor también tiene unidades de longitud como los lugares donde se sitúan los estratopuertos o el punto de inspección. Es por ello, que es conveniente que sean coherentes las unidades entre ellos. También es necesario que este valor sea mayor que 0. Aunque en este caso, es el algoritmo de las curvas de Dubins quien se asegura de comprobar que esta condición se cumpla. • Diferencial de tiempo o timestamp. Esta información corresponde con el paso de tiempo que el usuario quiera tener entre los puntos. Como para una velocidad constante y fijada no se puede tener un stepsize y un timestamp concretos, se van a calcular primero las trayectorias utilizando el stepsize correspondiente para después ser modificadas por el algoritmo para que respeten el timestamp indicado. De esta manera se tendrá un vector de tiempo más controlado aunque las distancias entre los puntos no sea siempre la misma. Esta situación es preferible ya que los vehículos no tripulados suelen tener una frecuencia de trabajo a la que reciben los datos. Por tanto, se le puede atribuir a este dato de tiempo el valor correspondiente a esta frecuencia. Por lo general, y como es lógico, este valor es siempre mayor que 0. Esta condición se comprueba en el algoritmo, dando un mensaje de error en el caso de que no se cumpla. • Visualización. En este argumento de entrada tiene que introducirse una variable tipo booleano, es decir, que tome valores true/f alse . En caso afirmativo o true, se haría una representación de las trayectorias que se han calculado con una leyenda que indique los distintos elementos que se pueden encontrar en la figura. Por otro lado, si se le da un valor de false, no se obtendría 2.2 Caso de uso 1: Inspección de un punto 35 ninguna figura al ejecutar el algoritmo pero sí se seguirían obteniendo los resultados en forma de matriz de puntos. • Simulación. Al igual que el argumento anterior, esta variable también es de tipo booleano. En este caso, realizaría una simulación encima de la representación en la que se observaría a una aeronave recorrer la trayectoria dibujada. De esta forma podemos comprobar que la trayectoria tiene los objetivos que se buscaban. Además, la simulación contiene una representación del Field of View de la aeronave. Cabe destacar que en el caso de que la visualización no se haya realizado porque no se haya solicitado, tampoco se haría la simulación. Una vez se tienen las entradas de la función, se reserva la primera parte del código para guardar estos datos en nuevas variables para tratar con ellos. Además, se crea el argumento de salida que corresponderá con una estructura denominada results . Un arreglo de estructura es un tipo de dato que agrupa datos relacionados mediante contenedores de datos llamados campos. Cada campo puede contener cualquier tipo de datos. Para acceder a los datos de una estructura, es posible usar una notación punto con el formato structName.fieldName [7]. Cálculo de la trayectoria mínima Tras el análisis de los datos iniciales introducidos, el algoritmo prosigue realizando el cálculo de las trayectorias con el formato conocido en 2 dimensiones de [x,y,θ] . Donde como siempre, x corresponde con la coordenada en el eje x e y es la coordenada sobre el eje y. Por último, θ es la orientación de la aeronave en el punto concreto. Para hacer este cálculo, primero hay que hallar cual es el camino más corto que puede realizar una aeronave sin ninguna especificación de tiempo. Esto se consigue con lo expuesto en el Código 2.11. Código 2.11 Cálculo de la trayectoria más corta. 1% Minimum path 2s = dl + 0.001; 3pAux = [pIns(1) - s*cos(pIns(3)), pIns(2) - s*sin(pIns(3)), pIns(3)]; 4pAux_FoV_post = [pIns(1) - dt*cos(pIns(3)), pIns(2) - dt*sin(pIns(3)), pIns(3)]; 5pAux_FoV_pre = [pIns(1) - dl*cos(pIns(3)), pIns(2) - dl*sin(pIns(3)), pIns(3)]; 6if dt > 0 7p = [pAux_FoV_pre; pAux_FoV_post; pIns; pAux; pAux_FoV_pre]; 8elseif dt < 0 9p = [pAux_FoV_pre; pIns; pAux_FoV_post; pAux; pAux_FoV_pre]; 10 else 11 p = [pAux_FoV_pre; pIns; pAux; pAux_FoV_pre]; 12 end 13 PATH = []; 14 for j = 1:length(p)-1 15 path = dubins_curve(p(j,:), p(j+1,:), r, stepsize, true); 16 PATH = [PATH; path]; 17 end La variable s que se encuentra en la línea 2 corresponde con la longitud del tramo de aproximación del circuito de espera que realiza la aeronave. Esta variable se modificará posteriormente para controlar el tiempo que tarda la aeronave en dar una vuelta. Para entender esta variable hay que entender los diferentes puntos auxiliares o waypoints que se usarán en el cálculo de la trayectoria. 36 Capítulo 2. Algoritmo Figura 2.28 Ejemplos puntos auxiliares (1). En la Figura 2.28 se observa un circuito de espera genérico. Los símbolos corresponden con los puntos auxiliares y el punto de inspección. El cuadrado naranja central es el punto de inspección, el triángulo amarillo de la izquierda es el punto denominado pAux_FoV_pre y el triángulo morado de la derecha es el punto denominado pAux_FoV_post . El triángulo amarillo indica en qué momento el punto de inspección entra en el rango de visión de los sensores. Por otro lado, el triángulo morado indica en que momento el punto de inspección sale del campo de visión. Esto nos indica que el valor de dt es negativo y, por tanto, el campo de visión incluye zonas detrás del paso de la aeronave. Lo que lleva a que se tenga que sobrepasar el punto de inspección para realizar una visualización completa del objetivo durante todo el tiempo posible, de ahí la situación del punto morado. Una vez quede esto claro, se puede entender la estructura condicional de las líneas 6-12. Esta líneas que se proponen, comparan si la longitud al borde posterior dt es mayor, menor o igual a 0. De esta manera sabemos si cuando el punto de inspección sale del rango de inspección, nos encontramos por detrás, por delante o sobre el mismo punto a inspeccionar respectivamente. Por tanto, sabiendo como se sitúan los puntos según el campo de visión de la aeronave, los puntos se reorganizan para pasar por todos ellos en orden. Por ejemplo, supóngase que dt es negativo y, por tanto, se encuentra pasado el punto de inspección, como en la Figura 2.28. Es por ello que el orden de los puntos sería el siguiente. Primero el punto correspondiente a la entrada del objetivo en el campo de visión, es decir, pAux_FoV_pre . En segundo lugar se pasaría por el punto de inspección. Y, en último lugar, se pasaría por el punto correspondiente a la salida del objetivo del campo de visión, es decir, pAux_FoV _post. Todavía no se ha mencionado ni explicado que posición corresponde al punto p_Aux . Este último será el controlado o desplazado por la variable s y lo hará de la siguiente manera. Si se observa cómo se han ordenado los vectores, todos ellos acaban pasando por el p_Aux y posteriormente por pAux_FoV_pre . Esto se debe a que estos son los puntos que en todos los casos se van a encontrar al inicio del circuito de espera y, por tanto, siempre hay que volver a ellos. En la línea 2 se define el valor de s como un poco más grande que dl , el cual es el la distancia al borde frontal. Es por ello que en el dibujo, Figura 2.28, el punto p_Aux se encontraría justo antes del triángulo amarillo, siendo este el nuevo punto más a la izquierda de trayectoria de aproximación. Esto permite que p_Aux pueda alargarse a su antojo (en este caso no podría acortarse porque es la trayectoria con mínima distancia) sin sobrepasar los límites que marcan los otros dos puntos auxiliares, ya que se alarga en la dirección opuesta. Finalmente, de la línea 13 a la 17, se realiza un simple bucle donde se realizan curvas de Dubins entre los waypoints creados y guardados en el vector p . Los resultados de estos caminos de Dubins 2.2 Caso de uso 1: Inspección de un punto 37 se guardan en la variable PATH en formato [x,y,θ]. El Código 2.12 continúa con los cálculos de la trayectoria más corta posible según los parámetros dados. Este calculan el tiempo que se tarda en realizar la ruta a la velocidad propuesta por el usuario. Código 2.12 Tiempos de la trayectoria más corta. 18 % Parameters minimum path 19 dist_min = length(PATH) * stepsize; 20 t_min = dist_min/v; 21 fprintf(’Minimum task time in seconds is %f \n’, t_min); 22 23 % Usuario check 24 prompt = ’Would you like to keep on with the run? (Yes = 1 | No = 0)\n’; 25 x = input(prompt); 26 if ~(x == 0 || x == 1) 27 error(’Error. Answer must be 0 or 1’) 28 elseif x == 0 29 return 30 end Para obtener el tiempo que se tarda en recorrer el camino más corto según los parámetros introducidos, primero hay que calcular la longitud de la trayectoria. Esta longitud se calcula de manera sencilla, ya que se conocen el número de puntos que componen la trayectoria y la distancia entre ellos, el stepsize. Multiplicando ambos valores, se obtiene la longitud de la trayectoria mínima (línea 2). Posteriormente, para obtener el tiempo que se tarda en recorrer esta distancia a una velocidad constante, utilizamos la ecuación v=e·t . En esta ecuación, v representa la velocidad del objeto, e es la distancia que se recorre y t el tiempo que se tarda en recorrer esa distancia. Se despeja el tiempo y se calcula su valor (línea 3). Las siguientes líneas de código corresponden a una verificación para que el usuario elija si quiere cambiar algún dato. Primero se informa por pantalla de cuál será el tiempo que se tarda en recorrer el trayecto mínimo. Más adelante, se le da a elegir al usuario si quiere seguir con la ejecución del algoritmo. Esto se ha introducido ya que una vez se sabe el tiempo mínimo que se tarda en dar una vuelta, el usuario tiene una idea de si le interesa introducir un tiempo superior o inferior. Este último haría que tuviera que introducirse un UAV adicional. Cálculo del número de UAVs necesarios y sus zonas de despegue A continuación se muestra el código correspondiente al cálculo del número de aeronaves necesarias para cumplir las especificaciones, Código 2.13. Estas especificaciones son la posición de los estratopuertos, la posición del punto de inspección, la velocidad de la aeronave, el radio de giro de la misma y la periodicidad de las inspecciones. Todos estos parámetros influyen a la hora de identificar cuántos UAVs son necesarios. Código 2.13 Cálculo del número de UAVs necesarios. 31 % Inspection frequency constraint 32 if t_freq ~= 0 33 dist_new = v*t_freq; 34 s = (dist_new - dist_min)/2; 38 Capítulo 2. Algoritmo 35 end 36 37 % s<0? Need more UAVs 38 nUAV = 1; 39 t_aux = t_freq; 40 while s < 0 41 t_aux = (nUAV + 1) * t_aux; 42 dist_new = v*t_aux; 43 s = (dist_new - dist_min)/2; 44 nUAV = nUAV + 1; 45 end Lo primero que se hace en esta parte es calcular cuál es la distancia que se recorrería si se volara durante un tiempo igual a la frecuencia de inspección. A esta distancia se la ha denominado dist_new (línea 3). En la siguiente línea de código se observa una variación del valor del parámetros s . El nuevo parámetro s corresponde con la diferencia de distancias entre la nueva variable dist_new y la dist_min calculada anteriormente. Esta resta corresponde con la distancia que falta por recorrer para alcanzar el tiempo de vuelta por inspección. Este valor se divide entre 2 porque si alargamos el tramo de aproximación, también se alarga el tramo contrario del circuito. Figura 2.29 Ejemplos puntos auxiliares (2). Lo comentado anteriormente se representa en la Figura 2.29. Los triángulos negros corresponden con los puntos ya conocidos pAux_FoV_pre y pAux_FoV_post . El triángulo amarillo es el nuevo punto auxiliar donde s es un valor mayor al mínimo visto en el Código 2.11. Como se puede observar, en rojo queda representada la distancia que aumenta la trayectoria respecto a la trayectoria mínima. Esta distancia aumenta el doble de lo que aumenta s porque afecta al tramo superior y al inferior del circuito. En el caso en el que el nuevo parámetro s sea mayor que 0, no hay ningún problema. Ya que, el tramo se puede alargar tanto como quiera para conseguir respetar la condición de tiempo. La dificultad aparece cuando el valor nuevo de s es menor que 0, es decir, la trayectoria mínima es más larga de lo que debería de ser si se quisiera cumplir la restricción de tiempos. La solución propuesta para esta cuestión es la que se muestra en las líneas 7-15 del Código 2.13. Dónde en el caso en el que s sea negativo, se aumenta en uno el número de UAVs necesarios y, por tanto, la trayectoria mínima se recorre en la mitad de tiempo. Más adelante se mostrará como se introducen los UAVs para que la separación entre ellos sea una fracción de la trayectoria total. Una 2.2 Caso de uso 1: Inspección de un punto 39 vez aumentado el número de UAVs se vuelve a hacer el cálculo del parámetros s y si se comprueba si ya es positivo o hay que añadir más aeronaves. Una vez conocido el número de aeronaves que van a ser necesarias para ejecutar la trayectoria respetando todas las restricciones impuestas, se procede a analizar si hay suficientes estratopuertos desde los que despegar, pudiéndose usar un mismo estratopuerto para más de un UAV. Aun así, en caso de que no haya spots suficientes para despegar, el algoritmo se encarga de repetir los estratopuertos que ya tiene registrados para que se utilicen en el despegue de más de una aeronave. Esto se puede observar a continuación en el Código 2.14. Código 2.14 Reorganización de zonas de despegue. 46 % UAVs take off spots 47 fprintf(’You need %d HAPS to achive the time constraint\n’, nUAV); 48 if size(p1,1) < nUAV 49 fprintf(’Not enough take off spots\n’) 50 prompt = ’Would you like to use the same spot for each UAV? (Yes = 1 | No = 0)\n’; 51 x = input(prompt); 52 if ~(x == 0 || x == 1) 53 error(’Error. Answer must be 0 or 1’) 54 elseif x == 0 55 return 56 end 57 while size(p1,1) < nUAV 58 p1 = [p1; repmat(p1,1)]; 59 end 60 end 61 62 % Order take off spots by closets to inspection 63 d = zeros(size(p1,1),1); 64 for i = 1:size(p1,1) 65 d(i) = dist(p1(i,:), pIns); 66 end 67 [~, I] = sort(d); 68 p1 = p1(I,:); Como se puede observar, el código le pregunta al usuario si quiere utilizar los estratopuertos indicados en las entradas de manera reiterada para todos los UAVs necesarios. Si el usuario responde afirmativamente, se repiten las coordenadas de las zonas de despegue tantas veces como UAVs sean necesarios. Por último en este trozo de código, se ordenan los estratopuertos por cercanía al punto de inspección, de esta manera se minimiza la distancia recorrida por los UAVs y se elije el estratopuerto más cercano. Cálculo de la trayectoria total en 2D A continuación, en el Código 2.15 se va a mostrar el cálculo de la trayectoria total. Código 2.15 Cálculo de la trayectoria 2D (1). 69 % Path calculation 70 pAux = [pIns(1) - s*cos(pIns(3)), pIns(2) - s*sin(pIns(3)), pIns(3)]; 40 Capítulo 2. Algoritmo 71 if dt > 0 72 p = [pAux_FoV_pre; pAux_FoV_post; pIns; pAux; pAux_FoV_pre]; 73 elseif dt < 0 74 p = [pAux_FoV_pre; pIns; pAux_FoV_post; pAux; pAux_FoV_pre]; 75 else 76 p = [pAux_FoV_pre; pIns; pAux; pAux_FoV_pre]; 77 end 78 n_laps = round(t_time/t_freq); % Inspection time constraint 79 fprintf(’UAV is going to do %d laps\n’, n_laps); 80 PATH = []; 81 trajectory_2D = struct; 82 for i = 1:n_laps 83 for j = 1:(length(p) - 1) 84 path = dubins_curve(p(j,:), p(j+1,:), r, stepsize, true); 85 PATH = [PATH; path]; 86 end 87 end Una vez se sepan el número de UAVs necesarios, las zonas de las que despegarán y los waypoints por los que pasarán, se comienza el desarrollo del cálculo de la trayectoria total. Este cálculo se va a realizar por tramos. Primero se va a calcular de manera independiente la trayectoria de inspección que realizará la aeronave. Para calcular esta trayectoria, se hará de la misma manera que el cálculo de la trayectoria de inspección mínima ya mostrada. A diferencia de que el punto auxiliar p_Aux , se encuentra en su nueva posición según el parámetro s . Posteriormente, se calculará el número de vueltas a dar según el tiempo por vuelta y el tiempo total de vuelo que introdujo el usuario. Finalmente, se creará un bucle donde se encadenen los waypoints creados en el orden correcto y se den tantas vueltas como se hayan calculado. Todos los resultados serán guardados en la variable PATH por el momento. Este PAT H será común para todas las aeronaves, pero el resto de la trayectoria no tiene por qué coincidir. Por tanto, a continuación en el Código 2.16 se realiza un bucle que calcule el resto de la trayectoria para cada UAV. Código 2.16 Cálculo de la trayectoria 2D (2). 88 for i = 1:nUAV 89 % Initial 90 path_ini = dubins_curve(p1(i,:), p(1,:), r, stepsize, true); 91 trajectory_2D(i).path_tot = path_ini; 92 trajectory_2D(i).path_length = length(path_ini); 93 trajectory_2D(i).path_time = length(path_ini)*stepsize/v; 94 95 % Holding pattern 96 trajectory_2D(i).path_tot = [trajectory_2D(i).path_tot; PATH]; 97 trajectory_2D(i).path_length = [trajectory_2D(i).path_length; length (PATH)]; 98 trajectory_2D(i).path_time = [trajectory_2D(i).path_time; length( PATH)*stepsize/v]; Como ya se ha mencionado, se realiza un bucle para cada UAV donde se guarden las trayectorias de cada uno. Las primeras trayectorias que se crean en el bucle, corresponden a los tramos del 2.2 Caso de uso 1: Inspección de un punto 41 camino total que van desde los estratopuertos hasta el primero punto dentro del circuito de espera. Este punto, denominado pAux_FoV_pre , es el que se encuentra en el instante donde el objetivo o punto de inspección, aparece por primera vez en el rango de visión. Esta trayectoria inicial se guarda en la estructura llamada tra jectory_2D en el campo path_tot . Se crean dos campos más llamados path_length y path_time donde se guardarán las longitudes de las trayectorias y el tiempo que se tarda en recorrer los tramos, respectivamente. Todas las trayectorias posteriores, se irán acumulando en esta estructura una tras otra. El siguiente tramo a guardar corresponde con la trayectoria ya calculada del circuito de espera. Las trayectorias se guardan de la misma manera que se han guardado las iniciales. Guardando además de la trayectoria en sí, su longitud y el tiempo que se tarda en recorrer. Se sabe que la trayectoria en el circuito de espera, acaba en el punto donde empezó, es decir, pAux_FoV_pre . Por tanto, es conveniente realizar una pasada más sobre el punto de inspección antes de retirarse de la misión. Por lo que se creó este tramo de inspección final con ese objetivo, Código 2.17. Código 2.17 Cálculo de la trayectoria 2D (3). 99 % Final Inspection 100 PATH_ins = []; 101 if dt > 0 102 for j = 1:2 103 path_ins = dubins_curve(p(j,:), p(j+1,:), r, stepsize, true); 104 PATH_ins = [PATH_ins; path_ins]; 105 end 106 elseif dt < 0 107 for j = 1:2 108 path_ins = dubins_curve(p(j,:), p(j+1,:), r, stepsize, true); 109 PATH_ins = [PATH_ins; path_ins]; 110 end 111 else 112 j = 1; 113 path_ins = dubins_curve(p(j,:), p(j+1,:), r, stepsize, true); 114 PATH_ins = [PATH_ins; path_ins]; 115 end 116 trajectory_2D(i).path_tot = [trajectory_2D(i).path_tot; PATH_ins]; 117 trajectory_2D(i).path_length = [trajectory_2D(i).path_length; length (PATH_ins)]; 118 trajectory_2D(i).path_time = [trajectory_2D(i).path_time; length( PATH_ins)*stepsize/v]; Al igual que en momentos anteriores, cuando se hace uso del vector que contiene los waypoints, es necesario verificar el orden de estos puntos auxiliares según el dato dt . Posteriormente, se realizan bucles que acaben cuando se haya revisado por completo el punto de inspección. A continuación de esto, se guarda la trayectoria de la misma manera que se ha hecho hasta ahora y se procede a comenzar el tramo de vuelta al punto de inicio o tramo final, Código 2.18. Código 2.18 Cálculo de la trayectoria 2D (4). 119 % Final 48 Capítulo 2. Algoritmo que estos casos de uso se han elegido porque tienen bastante aplicación en el ámbito de los UAVs o HAPS. Por ejemplo, se habló en la inspección del punto sobre la vigilancia de una zona pequeña como podría ser un campo de placas solares o un terreno agrario. Otras posibles aplicaciones que pueden tener estas trayectorias serían por ejemplo una misión de búsqueda y rescate o ampliar las comunicaciones en una zona más aislada. Póngase el primer caso como ejemplo. Un caso muy frecuentado en España son los incendios forestales. Estos ocurren sobre largas extensiones de tierra como puede ser la zona de Castilla y León, Extremadura o Aragón. Cuando se produce un incendio, cabe la posibilidad de que haya personas atrapadas en él o simplemente quiera observarse las dimensiones del mismo. Para ello, sería conveniente tener capacidad de monitorización de la zona y que las aeronaves puedan despegar desde distintos sitios, intentando elegir siempre el más cercano. Además, según el campo de visión de las aeronaves, estas tendrán que recorrer mayores distancias o menos. Por lo que también sería conveniente, distribuir toda la zona o área a revisar en distintas secciones las cuales sean vigiladas por distintos UAVs. Todas estas ideas, serán aplicadas a lo largo del algoritmo de inspección de un área. A continuación, y tal como se hizo en la sección anterior, se mostrará una parte del proceso y de los objetivos que se fueron marcando. Estos objetivos moldearon y le dieron forma a la trayectoria final. Además, como pasó con la inspección de punto, se impondrán una serie de constricciones que haya que respetar a la hora de realizar las rutas. La inspección del área, se planteó cuando el caso de un punto concreto estaba ya bastante desarrollado. Esto facilitó la velocidad de funcionamiento, ya que ya se tenía bastante experiencia trabajando con las curvas de Dubins. El primero de los objetivos o de las intenciones que se tenían con esta inspección de área, era recrear la misión que se ejecutaba en el artículo [19]. En esa misión, se comenzaba con un área circular que tenía que revisarse. Además, se disponía de un punto inicial desde el cual la aeronave despegaba. Una vez despegase, su propósito era acercarse a la zona por la parte más al Este en este caso y recorrer el círculo de inspección de Norte a Sur y viceversa alternando tramos rectilíneos. Para ello, las curvas de Dubins ayudaban bastante en la tarea, ya que el objetivo principal consistía en encontrar los waypoints de entrada y salida del área a inspeccionar. Los tramos paralelos se separaban entre sí, gracias al campo de visión de la aeronave. De esta manera, se aprovecha todo el Field of View posible, aunque existía un pequeño solapamientos entre pasadas controlado por un factor alpha. Todo esto se encuentra representado en la Figura 2.33. Figura 2.33 Ejemplo búsqueda en área (1). El método utilizado para resolver este problema, queda detallado en el artículo y consiste en 2.3 Caso de uso 2: Inspección de un área 49 calcular los puntos de entrada al área y de salida de ella en una pasada de Norte a Sur, o viceversa. Además, el algoritmo incluye una posible rotación de la trayectoria en el caso de que al usuario le interese, Figura 2.34. Figura 2.34 Ejemplo búsqueda en área (2). El código correspondiente a este cálculo se explicará cuando se comente el algoritmo. También, se probó la creación de estas trayectorias para mapas con diferentes formas. Ya que un mapa con forma de circunferencia no introducía muchas dificultades y no daba pie a hallar errores en el código. Por ejemplo, en la Figura 2.35 se muestran dos áreas distintas. La primera tiene forma de rectángulo y la segunda tiene la forma del contorno de Andalucía. Se puede comprobar que en ambas se ejecuta el código correctamente. (a) Área rectangular. (b) Área Andalucía. Figura 2.35 Otros ejemplos de búsqueda en área. Finalmente, cuando se consiguió reproducir y mejorar el procedimiento hallado en el artículo mencionado, se comenzó a pensar en distintas situaciones que pudieran enfocar el problema de otra manera según las necesidades del mismo. El siguiente objetivo, de la misma manera que sucedió en la inspección del punto, tiene que ver con cerrar la ruta. Ya que, cuando se termina de ejecutar la trayectoria, lo lógico es que el lugar final corresponda con un estratopuerto o una zona adecuada para el aterrizaje. Es por ello, que se plantearon distintas formas de cerrar la ruta. Sin embargo, utilizar el punto final que se tiene en la Figura 2.33 para desde ahí volver al inicio, no parece ser una solución muy óptima. Por tanto, la 50 Capítulo 2. Algoritmo solución propuesta fue modificar de manera sucinta los puntos de entrada y salida, para que dejaran un pasillo en la zona inferior del área sin vigilar. La aeronave aprovecharía la ruta de vuelta para recorrer esa zona sin inspeccionar hasta el momento, Figura 2.36. Figura 2.36 Ejemplo búsqueda en área cerrada. Una vez conseguido con éxito este propósito de cerrar la ruta y aprovechar el camino de vuelta. Se planteó como podría interesar introducir más de una aeronave, para los casos en los que el área sea tan grande que un solo UAV no consiga vigilar la zona con mucha rapidez. Una de las primeras ideas que surgieron, fue recorrer el camino total introduciendo un UAV tras otro, manteniendo una distancia de seguridad entre ellos. De esta manera, se podían dar vueltas siguiendo la trayectoria, y siempre estaría la mayor parte de la misma vigilada por aeronaves. Cuantas más aeronaves se tuviera, menos zonas no monitorizadas se tendría. Esta filosofía, es parecida a la que se empleó en el cálculo de las trayectorias para la inspección del punto. Aunque en esta ocasión, una vez conseguido el propósito, se observó que esa podría no ser la mejor solución al problema. Por tanto, se pensó en otra posible respuesta a la pregunta. Otra de las ideas fue comprobar primero la extensión del mapa, si este fuera demasiado grande en comparación con el campo de visión de la aeronave, se dividiría en secciones. Esta idea resultó ser mejor para vigilar más cantidad de terreno en menos tiempo. Ya que, el terreno total se divide en pequeñas porciones. Un inconveniente que se encontró fue dividir terrenos con formas irregulares. Esto desembocaba en una mayor dificultad y, por tanto, se optó por lo siguiente. En una primera estancia, el área a vigilar se sustituiría por un rectángulo que tenga la zona a inspeccionar inscrita en él. Para construir este rectángulo se haría uso de los puntos del área que se encuentren más al Norte, Sur, Este y Oeste. Una vez se tenga el rectángulo mencionado, se procede a dividir este rectángulo en secciones según las dimensiones del campo de visión del UAV. De esta manera, el usuario se asegura de lo largo que será el recorrido que haga cada UAV. Posteriormente, una vez se sepa cuantas divisiones de la zona se tienen, se manda a una aeronave a cada zona. Estas aeronaves saldrán de los estratopuertos más cercanos a su zona correspondiente y realizarán la inspección todas las aeronaves al mismo tiempo como se observaba en la Figura 2.36. De esta manera se aprovechará al máximo las trayectorias realizadas. A continuación se comentará el código que resuelve en MATLAB este problema y se irán detallando los objetivos descritos en este apartado. 2.3 Caso de uso 2: Inspección de un área 51 2.3.2 Desarrollo en MATLAB Argumentos de entrada Como se hizo en la inspección de punto, se irá explicando el código paso por paso indicando en qué consiste cada trozo del mismo y lo que hace. Primero se va a comenzar con los argumentos de entrada o inputs. La mayoría de ellos coinciden con las entradas del otro caso, pero se volverá a comentar su uso. • Estratopuertos. Como ya se mencionó en el otro algoritmo, los estratopuertos sirven como puntos de despegue y aterrizaje para las aeronaves que van a realizar esta misión. Se pueden introducir todos los que el usuario quiera y el algoritmo se encargará de organizarlos y seleccionar el más conveniente para cada caso. De la misma manera que en la inspección de punto, este dato se introduce como una matriz de 3 columnas con el formato [x,y,θ] . En cada fila se introducirá un nuevo estratopuerto. • Área de inspección. Este tipo de dato es nuevo en el trabajo, ya que en el caso anterior solo se proporcionaba un punto. En este caso, la entrada debe ser una matriz de 2 columnas que corresponde a las coordenadas en el eje x y en el eje y respectivamente. Esta matriz tiene que comenzar y acabar en el mismo punto para que se pueda considerar un área cerrada. Además, en los inicios del código, era necesario que esta matriz fuera bastante precisa y tuviera muchos puntos con poca distancia entre ellos. Pero con la evolución del algoritmo, eso ya no es necesario porque el propio código se encarga de crear un área adecuada para el análisis. • Velocidad de la aeronave. Mismo dato que en el caso anterior. Esta velocidad será utilizada a la hora de introducir la variable temporal. También, cabe destacar que no se admitirán valores negativos o 0 en el algoritmo. Por pantalla, se mostrará un mensaje de error en caso de que esto suceda. • Radio de giro mínimo de la aeronave. Otra vez se repite este dato que es básico para realizar los caminos de Dubins correspondientes a la trayectoria. No es necesario, como ya se comentó, que sea el radio mínimo de giro que pueda realizar la aeronave por su resistencia estructural o por temores de entrar en pérdida. También hay que comentar que, al igual que con la velocidad, no se permitirán valores que no sean positivos. • Tiempo total de vuelo. En este caso, las restricciones temporales difieren un poco de la inspección del punto. Esto es debido a que no se pretende recorrer el área entera en un cierto tiempo. Por tanto, esta variable indicará cual es el tiempo total de vuelo que quiere el usuario que esté la aeronave realizando su misión. Con este valor, se hará un cálculo del número de vueltas a realizar sobre el área. Como es lógico, un valor que no sea positivo de tiempo mostrará un mensaje de error diciendo que esto no es posible. • Campo de visión. El Field of View en esta ocasión, toma una gran importancia. Esto es debido a que, según las dimensiones de este, se realizará una división del área total en sub-áreas más grandes o más chicas. Además, también afecta a la separación lateral entre pasadas paralelas a lo largo del área. Estos conceptos se verán más adelante en la explicación de sus códigos correspondientes. Por último, en este caso se introduce un nuevo dato dentro de este vector. Hasta ahora, el vector incluía los valores del campo de visión: [dw,dt,dl] . En este caso, se va a utilizar también un dato que se va a denominar factor de solapamiento o alpha . Este dato se introducirá con los anteriores como una componente más del Field of View. Este factor de solapamiento servirá, como su propio nombre indica, para saber cuanta parte del campo de visión solapará 52 Capítulo 2. Algoritmo una pasada hecha con anterioridad. El valor de este factor debe encontrarse entre 0 y 1. Igual que en los ejemplos anteriores, en los datos del Field of View no se aceptarán valores negativos de la longitud hasta el borde superior (dl) o del ancho del campo de visión (dw). • Altitud. Esta constante, se utilizará de la misma manera que en el algoritmo de la inspección de punto. Se introducirá posteriormente como datos adicional a la trayectoria para que pueda realizarse una simulación más realista y algunas pruebas en un entorno real. Tampoco admitirá valores negativos, dando un mensaje de error en el caso en el que suceda. • Stepsize. El stepsize en esta ocasión tiene el mismo propósito que en los casos anteriores. Sirve como distancia que existirá entre punto y punto a la hora de crear las trayectorias en 2D. Aunque estos puntos serán sustituidos por otros que mantengan entre ellos un tiempo determinado conocido como timestamp. En este caso, si no se introduce un valor adecuado del mismo, el propio código de los caminos de Dubins se encargará de corregirlo. • Timestamp. Como se acaba de comentar, el timestamp es el tiempo correspondiente al cambio de las trayectorias al 4D. Este tiempo se utilizará convenientemente según la frecuencia a la que trabajen los vehículos no tripulados. Valores negativos no serán adecuados para esta variable temporal. • Visualización. Funcionando de la misma manera que se hizo en el caso anterior, esta variable de tipo booleano servirá para que el usuario decida si quiere una representación de la trayectoria sobre el mapa o no. Si el usuario quisiera una simulación, es obligatorio que marque esta variable como true para poder obtenerla. Ya que, no tiene sentido hacer una simulación sin una representación. • Simulación. Esta variable también es de tipo booleano como la anterior. En caso de que el usuario así lo decida, se realizará una simulación de las trayectorias recorridas por las aeronaves que sean necesarias. Estas serían todas las entradas correspondiente al nuevo algoritmo de la inspección del área. Una vez se hayan guardado estos argumentos en las variables correspondientes y haber modificado las unidades, por ejemplo de los ángulos, se procede a comenzar con el algoritmo de resolución del problema. Cálculo de área y sus divisiones En este apartado del código se va a convertir, como ya se comentó con anterioridad, el área a inspeccionar en un rectángulo que incluya de manera inscrita esta zona objetivo. Una vez se haga eso, según las dimensiones de este rectángulo, se dividirá el mismo en secciones rectangulares más pequeñas a las cuales se les asignará un único UAV. Para poder hacer esto, primero se va a comentar una función que se encuentra separada del código total, Código 2.25. Esta función es la encargada de, a partir de los puntos extremos de la zona a investigar, crear el rectángulo exterior a la misma. A esta función se la ha denominado rec_aIns.m . Código 2.25 Función rectángulo de inspección. 1function ap = rec_aIns(xW, xE, yN, yS, stepsize) 2V_length = floor((yN-yS)/stepsize)+1; 3H_length = floor((xE-xW)/stepsize)+1; 4ap = [xW*ones(V_length,1), (yS:stepsize:yN)’; 5(xW:stepsize:xE)’, yN*ones(H_length,1); 6xE*ones(V_length,1), (yN:-stepsize:yS)’; 7(xE:-stepsize:xW)’, yS*ones(H_length,1)]; 2.3 Caso de uso 2: Inspección de un área 53 8end Esta función necesita de los argumentos de entrada ya mencionados. Estos argumentos corresponden con los siguientes. Primero se solicita la coordenada x situada más al Oeste del área a inspeccionar, después se pide la misma coordenada pero más al Este. Posteriormente, del eje y se necesitan los coordenadas más al Norte y más al Sur. Además, también es una entrada un valor que corresponderá al stepsize, es decir, la distancia entre punto y punto del área a crear. Una vez se hayan leído los datos, la función calcula las distancias vertical y horizontal que formarán el rectángulo total. A estas distancias se les ha denominado V_length y H_length . Finalmente, se realiza un rectángulo cuyas dimensiones sean las calculadas y su posición sea teniendo los argumentos de entrada como extremos. La cantidad de puntos será tanta como el stepsize lo permita. Volviendo al código general de la inspección del punto, encontramos que comienza con las líneas que se observan en el Código 2.26 realizando un cálculo de los puntos exteriores al área. Estos serán usados posteriormente para crear el rectángulo adecuado. Código 2.26 División del área (1). 1ap = aIns; 2xE = max(ap(:,1)); 3xW = min(ap(:,1)); 4yN = max(ap(:,2)); 5yS = min(ap(:,2)); 6 7% Area division 8c1 = 3.5; % How many rep dw 9c2 = 5; % How many rep dl-dt 10 nH_len = floor((xE-xW)/(dw*c1)); 11 xW_vec = xW*ones(nH_len+1,1); 12 xE_vec = xE*ones(nH_len+1,1); 13 for i = 1:nH_len 14 xW_vec(i+1) = xW+i*(xE-xW)/nH_len; 15 xE_vec(i) = xW+i*(xE-xW)/nH_len; 16 end 17 nV_len = floor((yN-yS)/((dl-dt)*c2)); 18 yS_vec = yS*ones(nV_len+1,1); 19 yN_vec = yN*ones(nV_len+1,1); 20 for i = 1:nV_len 21 yS_vec(i+1) = yS+i*(yN-yS)/nV_len; 22 yN_vec(i) = yS+i*(yN-yS)/nV_len; 23 end Los puntos extremos a cualquier área dada, coinciden con los máximos y mínimos en cada eje. Después de calcular estos puntos, se comienza el proceso de dividir el rectángulo mayor en rectángulos más pequeños. Las constantes que aparecen creadas en las líneas 8 y 9, corresponden con unas constantes modificables cuyos valores tienen el siguiente significado. La primera de las constantes decide el número de pasadas que se quiere dar con un mismo UAV de Norte a Sur o viceversa. Lo que corresponde con el número de veces que se repite el ancho del campo de visión en la subdivisión y aporta la longitud horizontal del nuevo rectángulo. La segunda de las constantes es el equivalente vertical del valor anterior. Es decir, es el número de veces que se 54 Capítulo 2. Algoritmo repite el largo del campo de visión. Este valor se calcula con la resta entre la longitud hasta el borde superior y la longitud hasta el inferior. Una vez se haya decidido el valor de estas constantes, se calculan las coordenadas de los puntos más al Oeste y al Este de estas nuevas subdivisiones, y se guardan todas en los vectores xW_vec y xE_vec . Se realiza la misma operación para obtener las coordenadas de los puntos más al Norte y al Sur de las nuevas subdivisiones, y estos se guardan en los vectores yN_vec eyS_vec. Código 2.27 División del área (2). 24 trajectory_2D = struct; 25 k = nH_len*nV_len; 26 k_aux = 1; 27 for i = 1:nH_len 28 for j = 1:nV_len 29 trajectory_2D(k_aux).area = rec_aIns(xW_vec(i), xE_vec(i), yN_vec(j), yS_vec(j), stepsize); 30 trajectory_2D(k_aux).coord = [xW_vec(i), xE_vec(i), yS_vec(j), yN_vec(j)]; 31 k_aux = k_aux + 1; 32 end 33 end 34 nUAV = k; 35 fprintf([’Total area has been divided in %d parts. ’ ... 36 ’Each of them is going to be inspected by one UAV.\n’], nUAV); A continuación, se procede a calcular y guardar en una estructura estos nuevos rectángulos. El método elegido es ir creando rectángulos con la función ya explicada rec_aIns de izquierda a derecha y de abajo a arriba, Figura 2.37. Esto ha sido posible utilizando los vectores calculados y dos bucles concatenados recorriendo sus longitudes. Figura 2.37 Ejemplo área rectangular dividida. Todas las áreas y coordenadas de los extremos que se van calculando, se introducen en una estructura denominada tra jectory_2D . También se le asigna una posición a cada subdivisión controlada por la variable k_aux . Finalmente, se muestra un mensaje por pantalla que dice en cuántas partes se ha dividido el área total y se comunica que todas serán recorridas por un único UAV. 2.3 Caso de uso 2: Inspección de un área 55 El siguiente paso será identificar las distintas zonas de despegue introducidas como argumentos de entrada. Posteriormente, se analizará la necesidad de incluir más o repetir alguna. El código correspondiente a lo comentado coincide con el que se realizó para la inspección del punto y está escrito en el Código 2.28. Código 2.28 Reorganización de zonas de despegue. 37 % UAVs take off spots 38 if size(p1,1) < nUAV 39 fprintf(’Not enough take off spots\n’) 40 prompt = ’Would you like to use the same spot for each UAV? (Yes = 1 | No = 0)\n’; 41 x = input(prompt); 42 if ~(x == 0 || x == 1) 43 error(’Error. Answer must be 0 or 1’) 44 elseif x == 0 45 return 46 end 47 while size(p1,1) < nUAV 48 p1 = [p1; repmat(p1,1)]; 49 end 50 end Al igual que en la inspección de un punto en el mapa, en el caso de que no hubiera suficientes zonas de despegue, se le preguntaría al usuario si quiere reutilizar las que ya se introdujeron. En este caso, se seguirían introduciendo las coordenadas de los estratopuertos a reutilizar, para posteriormente ser reorganizadas y elegidas según la proximidad a la ruta. Cálculo de la trayectoria total 2D En esta sección, se calcularán las trayectorias que tenga que realizar cada HAPS según el área que se le asigne. De esta manera, se asegura que cada zona a inspeccionar tenga una trayectoria calculada. Para hacer esto primero hay que elegir una subdivisión y guardar sus coordenadas, Código 2.29. Código 2.29 Cálculo de trayectoria 2D (1). 51 % Path calculation 52 for k_aux = 1:nUAV 53 xW = trajectory_2D(k_aux).coord(1); 54 xE = trajectory_2D(k_aux).coord(2); 55 yS = trajectory_2D(k_aux).coord(3); 56 yN = trajectory_2D(k_aux).coord(4); 57 ap = trajectory_2D(k_aux).area; 58 aprox_factor = (-1)*log10(stepsize); 59 count = 0; 60 61 % Initial 62 pAux = [xW, yS, 0]; 63 64 % Select take off spots by closets to initial point 65 d = zeros(size(p1,1),1); 56 Capítulo 2. Algoritmo 66 for i = 1:size(p1,1) 67 d(i) = dist(p1(i,:), pAux); 68 end 69 [~, I] = sort(d); 70 p_ini = p1(I(1),:); 71 72 path_tot = []; 73 PATH = dubins_curve(p_ini, pAux, R, stepsize, true); 74 trajectory_2D(k_aux).length_ini = length(PATH); 75 path_tot = [path_tot; PATH]; Una vez creado el bucle y guardado las coordenadas de la zona correspondiente, se elije el punto más al suroeste como punto inicial de la trayectoria. Posteriormente, se selecciona de toda la lista de zonas de despegue, la más cercana a este punto inicial. El siguiente paso será calcular la trayectoria inicial que va desde el estratopuerto más cercano a la zona hasta el punto más al suroeste de la misma. El siguiente paso es comenzar el bucle que vaya calculando pasadas de Norte a Sur y viceversa hasta llegar al final de la subdivisión. Código 2.30 Cálculo de trayectoria 2D (2). 76 while 1 77 p_lap = []; 78 i = 0; 79 80 % Inspection 81 while alpha*dw*(i+1/2) < xE-xW 82 xi = round(xW + (i+1/2)*alpha*dw,aprox_factor); 83 ind = [find(round(ap(:,1),aprox_factor) == xi)]; 84 if ap(ind(1),2) < ap(ind(end),2) 85 yni = ap(ind(end),2); 86 ysi = ap(ind(1),2); 87 else 88 yni = ap(ind(1),2); 89 ysi = ap(ind(end),2); 90 end 91 p_lap = [p_lap; xi ysi pi]; El método para recorrer el área es el que se utilizó en los primeros intentos de estas trayectorias. A continuación se va a explicar este método cuyo código queda escrito en el Código 2.30. Primero, se elije una coordenada en el eje x que se encuentre a una distancia equivalente a la mitad del ancho del campo de visión. Una vez se calcule esta coordenada, se busca en el área los puntos correspondientes a ese valor. Se obtendrá un punto en la parte superior del rectángulo y otro en la parte interior. La Figura 2.38 muestra un dibujo explicativo de lo planteado. 2.3 Caso de uso 2: Inspección de un área 57 Figura 2.38 Selección de coordenadas de pasadas de inspección (1). Como se ha comentado, los puntos en gris servirán como referencias para calcular los puntos de entrada y de salida a la zona de inspección. Cabe destacar que el punto que se encuentra en el Sur del área, se guardará en una matriz de puntos denominada p_lap con la que se regresará a la posición inicial como se explicó en la sección de introducción de este algoritmo. Los puntos de entrada y de salida, se calculan como viene descrito en el Código 2.31. Código 2.31 Cálculo de trayectoria 2D (3). 92 if mod(i,2) == 1 93 pAuxSi = [xi, yni+dl, -pi/2]; 94 pAuxEi = [xi, ysi-(1-alpha*dw/(4*dt))*dt, -pi/2]; 95 else 96 pAuxSi = [xi, ysi+(1+dw/(4*dl))*dl, pi/2]; 97 pAuxEi = [xi, yni-dt, pi/2]; 98 end 99 p = [pAux; pAuxSi; pAuxEi]; 100 for j = 1:length(p)-1 101 PATH = dubins_curve(p(j,:), p(j+1,:), R, stepsize, true); 102 path_tot = [path_tot; PATH]; 103 end 104 pAux = p(3,:); 105 i = i + 1; 106 end A los puntos de entrada a la zona se les ha denominado pAuxSi y pAuxEi , los nombres provienen de inicio, Start y fin, End. Los tramos se realizarán comenzando en el último punto calculado pAux , para posteriormente dirigirse al punto de inicio y acabar en el punto final. El punto final se convertirá tras terminar el bucle en el nuevo punto auxiliar pAux. El cálculo de los puntos de entrada y salida se realiza teniendo en cuenta el campo de visión. El objetivo es que se recorra toda la zona de inspección con un tramo rectilíneo de punta a punta encontrándose esta siempre en el rango de visión de la aeronave, al igual que se hizo con la inspección del punto. Una vez se tengan estos puntos, se construye una cadena que contenga las posiciones calculadas y se realizan curvas de Dubins entre los distintos puntos al igual que como se ha hecho hasta el momento. Este bucle de operaciones continuará hasta que se llegue al final de la zona a investigar, Figura 2.39. 64 Capítulo 2. Algoritmo De esta manera finaliza el segundo algoritmo correspondiente con las inspección del área. A continuación se explicará de forma breve como se han unido estos dos algoritmos para crear una sola función que identifique cuál de los dos es más conveniente para la situación que propone el usuario que introduce los datos. 2.4 Algoritmo completo Por último, en esta sección se van a comentar los distintos modos de uso de estas funciones. Estos modos se dividirán en ejecutar los algoritmos a través de scripts de MATLAB o hacer uso de una aplicación creada con MATLAB App Designer. Aunque antes de esto, será necesario unificar ambos algoritmos en una única función que seleccione cual de los dos va a ejecutar. Esta nueva función general recibirá como entradas una serie de datos guardados en una estructura. Estos valores corresponden a todas las variables necesarias para identificar cual sería la mejor ruta. Esta estructura se descompondrá guardando los valores en variables concretas que se introducirán como inputs en los propios algoritmos. Además, se creará una función dentro de la propia función general con el propósito de identificar cuál de los dos algoritmos es el que hay que utilizar. Para ello, se analizará el dato de la estructura correspondiente a la inspección. Según si este dato contiene una serie de puntos los cuales formen un área o si es un único punto a inspeccionar, se seleccionará un caso de uso u otro. 2.4.1 Función unificadora Ambos algoritmos de planificación de trayectorias se pueden ejecutar de manera individual introduciendo los datos solicitados en ellas. Hay que tener en cuenta que la entrada para la inspección de punto debe ser un único punto con las coordenadas correspondientes y que los datos introducidos de velocidad, tiempo, posición tengan unidades adecuadas y lógicas. Aunque puedan ejecutarse las funciones introduciendo los datos directamente en la ventana de comandos, es más conveniente realizar un archivo a parte donde se puedan modificar los datos e introducirlos en las funciones. Es por ello que se creó una nueva función de MATLAB. Esta función se la denominó inspection.m. Esta función recibe los datos que el usuario quiera introducir en forma de estructura de datos. Posteriormente, analiza estos datos para comprobar si los valores de las variables tienen las unidades correctas y son magnitudes lógicas. Finalmente, elige la función que más se adecua a los datos proporcionados y ejecuta el algoritmo correspondiente. Si el usuario lo solicitase, sacaría por pantalla una representación de las trayectorias, una simulación y guardaría los resultados en formato .csv para realizar pruebas en un entorno real. El Código 2.41 muestra la función que se ha comentado. Código 2.41 Algoritmo de inspección genérico (1). 1%%INSPECTION 2% This function will choose between the area searching task or the point 3% inspect task, basing on the initial data. Also it will save in csv 4% format the obtained trajectories. 5 6function results = inspection(data) 7% Select task 8selector = task_selector(data); 9if selector 2.4 Algoritmo completo 65 En esta parte del código se observa una breve descripción de lo que va a hacer la función además del inicio de la misma. En el comienzo de esta se puede observar que se hace uso de una función interna del algoritmo llamada task_selector . Esta función elige entre los algoritmos de la inspección de un punto o la inspección de un área cuál será ejecutado. La función selectora se analizará al final de esta sección. Código 2.42 Algoritmo de inspección genérico (2). 10 % Save Inputs 11 p1 = data.InitialPoints; 12 pIns = data.Inspectdata; 13 v = data.UAVspeed; 14 r = data.UAVradius; 15 t_freq = data.Taskfrequency; 16 t_time = data.Tasktime; 17 dw = data.FoV(1); 18 dt = data.FoV(2); 19 dl = data.FoV(3); 20 altitude = data.altitude; 21 stepsize = data.stepsize; 22 dtime = data.dtime; 23 doPlot = data.Plot; 24 doSim = data.Simulation; A continuación, Código 2.42, se guardan los datos de la estructura, introducida como único argumento de la función, en una serie de variables para analizar sus valores y emitir mensajes de error en caso de que sea necesario. Código 2.43 Algoritmo de inspección genérico (3). 25 % Check inputs 26 if ~(size(p1,2) == 3) 27 error(’Error. Initial points have non valid format’) 28 end 29 if ~(size(pIns,2) == 3) 30 error(’Error. Inspection point has non valid format’) 31 end 32 if ~(v > 0) 33 error(’Error. UAV speed must be positive’) 34 end 35 if ~(r > 0) 36 error(’Error. UAV turning radius must be positive’) 37 end 38 if ~(t_freq >= 0) 39 error(’Error. Task inspection frequency has non valid value’) 40 end 41 if ~(t_time > 0) 42 error(’Error. Task inspection time has non valid value’) 43 end 44 if t_time < t_freq 66 Capítulo 2. Algoritmo 45 error(’Error. Task inspection time must be greater than task inspection frequency’) 46 end Código 2.44 Algoritmo de inspección genérico (4). 47 if ~(dw > 0) 48 error(’Error. FoV width must be positive’) 49 end 50 if ~(dl > 0) 51 error(’Error. FoV leading edge must be positive’) 52 end 53 if ~(dl > dt) 54 error(’Error. FoV leading edge (dl) must be greater than FoV trailing edge (dt)’) 55 end 56 if ~(altitude > 0) 57 error(’Error. UAV altitude must be positive’) 58 end 59 if ~(stepsize > 0) 60 error(’Error. Stepsize must be positive’) 61 end 62 if ~(dtime > 0) 63 error(’Error. Timestamp must be positive’) 64 end 65 if ~(doPlot == true || doPlot == false || doSim == true || doSim == false) 66 error(’Error. Plot and Simulaton data must be logical value’) 67 end 68 if doPlot == false && doSim == true 69 error(’Error. Cannot do a simulation without the graphic representation’) 70 end Como se explicó en los algoritmos de ambos casos de uso, se comprueba que las entradas tengan valores lógicos, como son números positivos para variables de tiempo o velocidad, o vectores de tres coordenadas para indicar posiciones, Código 2.43. También se comprueba que las variables que se utilizan para indicar si se quiere una representación o una simulación sean variables de tipo booleano, Código 2.44. Código 2.45 Algoritmo de inspección genérico (5). 71 % Run task 72 fprintf(’The Point Inspect Task will be performed\n’) 73 results = point_inspect(p1, pIns, v, r, t_time, t_freq, data.FoV, altitude, stepsize, dtime, doPlot, doSim); Finalmente, tras comprobar que todos los valores de las entradas tienen valores adecuados, se procede a informar del algoritmo que se va a utilizar y se ejecuta el mismo, Código 2.45. De forma 2.4 Algoritmo completo 67 análoga, se hace lo mismo en el caso de la inspección del área, Códigos 2.46, 2.47, 2.48 y 2.49. A diferencia que en la inspección del punto, en el área no se comprueba ningún tiempo de vuelta ya que no se aporta como datos. Aunque si hay que verificar otras variables como el factor de solapamiento alpha. Código 2.46 Algoritmo de inspección genérico (6). 74 else 75 % Save Inputs 76 p1 = data.InitialPoints; 77 aIns = data.Inspectdata; 78 v = data.UAVspeed; 79 r = data.UAVradius; 80 t_time = data.Tasktime; 81 dw = data.FoV(1); 82 dt = data.FoV(2); 83 dl = data.FoV(3); 84 altitude = data.altitude; 85 alpha = data.FoV(4); 86 stepsize = data.stepsize; 87 dtime = data.dtime; 88 doPlot = data.Plot; 89 doSim = data.Simulation; Código 2.47 Algoritmo de inspección genérico (7). 90 % Check inputs 91 if ~(size(p1,2) == 3) 92 error(’Error. Initial points have non valid format’) 93 end 94 if ~(size(aIns,2) == 2) 95 error(’Error. Inspection area has non valid format’) 96 end 97 if ~(v > 0) 98 error(’Error. UAV speed must be positive’) 99 end 100 if ~(r > 0) 101 error(’Error. UAV turning radius must be positive’) 102 end 103 if ~(t_time > 0) 104 error(’Error. Task inspection time must be positive’) 105 end 106 if ~(dw > 0) 107 error(’Error. FoV width must be positive’) 108 end 109 if ~(dl > 0) 110 error(’Error. FoV leading edge must be positive’) 111 end 112 if ~(dl > dt) 68 Capítulo 2. Algoritmo 113 error(’Error. FoV leading edge (dl) must be greater than FoV trailing edge (dt)’) 114 end Código 2.48 Algoritmo de inspección genérico (8). 115 if ~(altitude > 0) 116 error(’Error. UAV altitude must be positive’) 117 end 118 if ~(alpha > 0 && alpha <= 1) 119 error(’Error. Lane spacing must be between 0 and 1’) 120 end 121 if ~(stepsize > 0) 122 error(’Error. Stepsize must be positive’) 123 end 124 if ~(dtime > 0) 125 error(’Error. Timestamp must be positive’) 126 end 127 if ~(doPlot == true || doPlot == false || doSim == true || doSim == false) 128 error(’Error. Plot and Simulaton data must be logical value’) 129 end 130 if doPlot == false && doSim == true 131 error(’Error. Cannot do a simulation without the graphic representation’) 132 end Código 2.49 Algoritmo de inspección genérico (9). 133 % Save trajectories as .csv 134 doCSV = data.SaveCSV; 135 if doCSV 136 for i = 1:length(results) 137 id = num2str(i); 138 name = ’trajectory_uav_XXX.csv’; 139 switch length(id) 140 case 1 141 name(16:17)=’0’; 142 name(18)=id; 143 case 2 144 name(16)=’0’; 145 name(17:18)=id; 146 case 3 147 name(16:18)=id; 148 end 149 csvwrite(name,results(i).trajectory); 150 end 151 end 152 end 2.4 Algoritmo completo 69 A continuación, se van a indicar las líneas de código correspondientes a la selección del algoritmo a utilizar, Código 2.50. Esto se realiza comprobando si el dato correspondiente al objeto a inspeccionar es un único punto o es una serie de puntos que conforman un área completa. Según el resultado de esta comprobación se le da un valor de 0 o 1 a una variable que se introducirá como condición de la estructura condicional mostrada en los códigos anteriores. Código 2.50 Algoritmo de inspección genérico (10). 153 function selector = task_selector(data) 154 aux = data.Inspectdata; 155 if size(aux,1) == 1 156 selector = 1; 157 else 158 selector = 0; 159 end 160 end Una vez explicado el código correspondiente al algoritmo completo, se va a exponer las maneras dos maneras que se han pensado para introducir los datos en esta función. Como ya se comentó al principio de la sección, una manera consistirá en crear un script de MATLAB donde se introduzcan los datos en una estructura de forma directa para posteriormente ejecutar la función en el propio archivo. La otra de las maneras consistirá en crear una aplicación con el App Designer de MATLAB como se comentará más adelante. 2.4.2 Entradas a través de un script La primera de las opciones mencionadas consiste en crear de forma individual una estructura que contenga los campos que se solicitan en la función inspection.m . Esta es la manera más sencilla para realizar pruebas, obtener resultados y ha sido la utilizada durante la creación del algoritmo. Para facilitar cuales son estos datos se ha creado un archivo .m de MATLAB que comenta y explica estos campos. En el Código 2.51 se muestran los datos mencionados. Código 2.51 Script de datos iniciales (1). 1% UAV speed: data.UAVspeed [m/s] 2% UAV turning radius: data.UAVradius [m] 3% Take Off spot: data.InitialPoints [x y yaw] [m, o ¯] 4% Field of View: data.FoV [dw dt dl alpha] [m] 5% Altitude: data.altitude [m] 6% Inspection data: data.Inspectdata [x y yaw] Punto [1x3] or área [mx3] 7% Total flight time: data.Tasktime [s] 8% Inspection frequency: data.Taskfrequency [s] 9% Graphic representation: data.Plot [true or false] 10 % Simulation of the results: data.Simulation [true or false] 11 % Save results as CSV: data.SaveCSV [true or false] 12 % Distance between points: data.stepsize [m] 13 % Time between WP: data.dtime [s] 14 15 clc, clear all, close all 70 Capítulo 2. Algoritmo Primero se crean unas líneas comentadas donde se explica cuales son los datos que hay que introducir, después se indica el formato para introducirlos en una estructura llamada data y, por último, se indican las unidades de cada entrada. Después de esto, se elimina por seguridad todos los datos, gráficas y líneas de código existentes para que no exista acumulación de datos antiguos y se produzcan errores. Código 2.52 Script de datos iniciales (2). 16 data.UAVspeed = 0.75; 17 data.UAVradius = 0.15; 18 data.InitialPoints = [0 -3.5 -90; -4.5 -2 0]; 19 data.FoV = [0.5 -0.3 0.3 0.9]; 20 data.altitude = 1; 21 data.Inspectdata = rec_aIns(-3, 3, 2, -2.5, 0.1); 22 data.Tasktime = 40; 23 data.Taskfrequency = 6; 24 data.Plot = true; 25 data.Simulation = false; 26 data.SaveCSV = false; 27 data.stepsize = 0.01; 28 data.dtime = 0.1; 29 30 results = inspection(data); Una vez introducidos, como se indicaba en los comentarios, todos los datos en la estructura y en los campos correspondientes, se ejecuta la función inspection.mexplicada en el apartado anterior la cual realiza las trayectorias que corresponden con la situación introducida, Código 2.52. 2.4.3 Entradas a través de MATLAB App Designer La segunda manera que se utiliza para el almacenamiento de los datos y la ejecución del algoritmo sería utilizar el módulo de App Designer de MATLAB. App Designer permite crear aplicaciones profesionales de manera cómoda e intuitiva. Este método facilita el entendimiento de los datos introducidos sin saber sobre la estructura de datos necesaria para la ejecución del código. Está planteada para futuros usos del mismo dónde otros usuarios no tengan porqué saber sobre las entrañas del algoritmo [1]. En este proyecto no se va a mostrar el proceso de creación de la aplicación, aunque si se comentarán los elementos que se han introducido y como se utilizan para obtener los datos solicitados por el algoritmo. El nombre que se le ha dado a la aplicación es InspectionApp.mlapp . Cuando este archivo se ejecuta, aparece una ventana como la que se muestra en la Figura 2.42. 2.4 Algoritmo completo 71 Figura 2.42 Aplicación del algoritmo completo (1). En la imagen se puede observar los distintos datos que se mencionaron en los puntos anteriores. En esta ocasión están expuestos como campos a los que introducirles los valores numéricos que el usuario quiera. Además, se añade un dibujo que muestran los parámetros correspondientes al campo de visión. Esto se ha decidido así para facilitar el entendimientos de los valores a introducir. Se pueden observar los campos correspondientes a la velocidad de la aeronave, el radio de giro, los parámetros que determinan el campo de visión y los parámetros de tiempos. También aparecen una serie de casillas donde se puede seleccionar si se quiere la representación de las trayectorias, la simulación de una aeronave sobre ellas y si se quiere obtener las trayectorias en formato .csv . Otros parámetros importantes son el stepsize y el timestamp, estos se encuentran a parte por su importancia y delicadeza a la hora de ser introducidos. Por último, se puede observar un cuadro y unas casillas para introducir coordenadas de estratopuertos. El método de utilización de esta sección consiste en añadir unas coordenadas con el formato indicado y pulsar el botón correspondiente. De esta manera se introduce el estratopuerto indicado de manera sencilla e intuitiva como se observa en la Figura 2.43. Figura 2.43 Aplicación del algoritmo completo (2). El módulo de App Designer crea un código correspondiente a todos los elementos que se van introduciendo en tu ventana de trabajo. Aunque también se pueden introducir líneas de código en 72 Capítulo 2. Algoritmo caso necesario. Por ejemplo, en esta ocasión ha sido necesario escribir unas cuantas líneas para asignar los valores que se introducen en la ventana inicial a las variables correspondientes en la estructura de datos, Código 2.43. También ha sido necesario añadir una función que ejecute el algoritmo completo cuando se pulse el botón "play". Y, por último, otra función fue necesaria para realizar la adición de los distintos estratopuertos, Código 2.54. Código 2.53 Código App Designer (1). 1% Button pushed function: PlayButton 2function PlayButtonPushed(app, event) 3 4% Get the values of each data 5data.UAVspeed = app.SpeedEditField.Value; 6data.UAVradius = app.TurningratiusEditField.Value; 7data.FoV = [app.WidthdwEditField.Value app.TrailingedgedtEditField. Value app.LeadingedgedlEditField.Value app.AlphafactorEditField. Value]; 8data.altitude = app.AltitudeEditField.Value; 9data.Tasktime = app.FlighttimeEditField.Value; 10 data.Taskfrequency = app.InspectiontimeEditField.Value; 11 data.Plot = app.PlotCheckBox.Value; 12 data.Simulation = app.SimulationCheckBox.Value; 13 data.SaveCSV = app.SaveCSVCheckBox.Value; 14 data.stepsize = app.StepsizeEditField.Value; 15 data.dtime = app.TimestampEditField.Value; 16 17 % Get initial points values 18 cell = app.SaveCoord.Value; 19 data.InitialPoints = str2double(strsplit(cell{2,1},’,’)); 20 for i = 3:length(cell) 21 data.InitialPoints = [data.InitialPoints; str2double(strsplit( cell{i,1},’,’))]; 22 end 23 % data.Inspectdata = rec_aIns(-3, 3, 2, -2.5, 0.1); 24 data.Inspectdata = [1 1 0]; 25 26 % Do inspection function 27 inspection(data); 28 end Código 2.54 Código App Designer (2). 1% Button pushed function: AddstratoportsButton 2function AddstratoportsButtonPushed(app, event) 3 4% Add new point to table 5newPoint = app.StratoportEditField.Value; 6text = "[" + newPoint + "]"; 7app.StratoportsCoord.Value = [app.StratoportsCoord.Value; text]; 2.4 Algoritmo completo 73 8app.SaveCoord.Value = [app.SaveCoord.Value; newPoint]; 9end Finalmente, cuando se hayan decido los datos iniciales, el usuario tiene que pulsar el botón con la palabra "play". Esto hará que se ejecute la función completa ya explicada con los datos introducidos. Los cuadros de diálogos correspondientes seguirán apareciendo en la ventana de comandos de MATLAB. Posteriormente se obtendrán las trayectorias y la representación correspondiente en caso de que se haya solicitado. Este formato en forma de aplicación es un prototipo inicial. Tiene muchas posibles mejoras en cuanto a su formato que se comentarán más adelante en las conclusiones. 80 Capítulo 3. Simulación 3.1.2 Ejemplo 2 En este segundo ejemplo de trayectoria para la inspección de un punto, se va a complicar un poco más la situación respecto al primer ejemplo. En esta ocasión, se introducirá un tiempo de vuelta menor al tiempo mínimo. Esto provocará que el algoritmo tenga que utilizar más de un UAV para cumplir esta condición de tiempo. Al igual que en el primer ejemplo, a continuación, en la Figura 3.6, se expondrán dos imágenes con los datos de la nueva trayectoria. La primera imagen corresponderá con el código correspondiente a estos datos y la segunda contendrá la estructura que se crea al ejecutar este código. (a) Código de los datos. (b) Workspace de los datos. Figura 3.6 Datos del ejemplo 2. En esta ocasión se han modificado los datos iniciales a excepción de las especificaciones de las aeronaves, que se han mantenido los datos de velocidad y radio de giro. El punto inicial es el mismo, pero el punto objetivo se ha modificado para que se encuentre en la posición [2,1] y tenga que ser inspeccionado con un ángulo de 30 grados respecto del eje x en sentido horario. El campo de visión y la altitud a la que volarán los vehículos también se han considerado los mismos. Las constricciones de tiempo se han modificado, ya que, como se observará en la ventana de comandos, se ha introducido un tiempo de vuelta menor que el tiempo mínimo para inspeccionar. Esto provocará que se tenga que añadir un segundo UAV. Además, el tiempo total de vuelo se ha introducido de manera que se den dos vueltas al circuito de inspección. En la Figura 3.7 se puede observar el mensaje que sale por pantalla cuando se ejecuta el algoritmo en la imagen de la izquierda. La gráfica de la derecha muestra la trayectoria final una vez contestadas las preguntas emergentes. (a) Mensajes emergentes. (b) Representación. Figura 3.7 Ruta completa del ejemplo 2. 3.1 Inspección de un punto 81 En este caso existen algunas diferencias respecto al ejemplo pasado. Se puede observar que de la misma manera que en el primer ejemplo, se va a realizar la inspección del punto. El tiempo mínimo de vuelta es superior al tiempo indicado, y por tanto, el algoritmo señala que se necesitarán dos aeronaves para poder cumplir la condición. Además, como en los datos únicamente se indicó un solo punto de despegue o estratopuerto, el algoritmo pregunta si se desea utilizar el mismo punto de despegue para ambos UAVs. Finalmente se indica que sí y se informa de que se darán dos vueltas. Si se observa la representación de la ruta que se va a realizar desde una vista aérea de la misma, se pueden ver los 30 grados de inclinación del circuito de espera respecto al eje x en sentido horario. Se ha omitido la leyenda para visualizar mejor toda la ruta, ya que los elementos son los mismos que en el ejemplo anterior. A continuación, se va a mostrar una serie de instantes, Figura 3.8y Figura 3.9, en los que se observará a la flota de dos aeronaves recorrer estas trayectorias con una distancia de seguridad entre ellas. (a) Despegue UAV 1. (b) Despegue UAV 2. (c) Inspección UAV 1. Figura 3.8 Despegues y aproximación. Ejemplo 2. (a) Inspección equidistanciada. (b) Regreso UAV 1. (c) Regreso UAV 2. Figura 3.9 Inspección y aterrizajes. Ejemplo 2. En este conjunto de figuras se encuentran seis momentos importantes a lo largo de la simulación. El primer momento correspondiente con el despegue de la primera aeronave mientras que la segunda permanece en espera. La segunda imagen muestra como el UAV 1 ya está aproximándose al punto de inspección mientras que el UAV 2 espera el momento exacto para realizar su despegue. Esta espera permitirá que ambas aeronaves mantengan una distancia de seguridad equivalente a la mitad de la longitud del circuito de espera. En la tercera imagen, el UAV 1 ya ha realizado la primera inspección y procede a recorrer el circuito de espera. El segundo UAV se encuentra en una aproximación al punto objetivo. La 82 Capítulo 3. Simulación cuarta imagen muestra a ambos UAVs recorriendo la ruta de inspección manteniendo una distancia específica como ya se ha comentado. Tras dar dos vueltas al circuito, el primer UAV regresa al estratopuerto para realizar su aterrizaje como se muestra en la quinta imagen. Finalmente, una vez aterrizado el primero UAV, se encuentra el UAV 2 de regreso al punto inicial. Al igual que en el primer ejemplo, una vez termina la simulación, el algoritmo guarda en una estructura todas estas trayectorias planteadas. En este caso se tendrán dos elementos dentro de la estructura que corresponderán a las dos aeronaves respectivamente. Además, se van a mostrar algunos de los datos de ambas aeronaves para comparar los resultados. Figura 3.10 Resultados obtenidos del ejemplo 2 (1). (a) UAV 1. (b) UAV 2. Figura 3.11 Resultados obtenidos del ejemplo 2 (2). Como ya se ha comentado, en la Figura 3.10 se observa la estructura que contiene las dos trayectorias. La fila 1 corresponde al primer UAV y la fila 2 corresponde al segundo. El UAV 2 tiene más elementos, como se puede comprobar. Esto es debido a que tiene que realizar una espera inicial mientras el primer UAV ya está recorriendo parte de la ruta. Esta situación que se menciona se puede observar en la Figura 3.11. La primera columna de las dos matrices de puntos corresponde con el tiempo de vuelo. En este caso las imágenes muestran las posiciones donde se encuentran ambos UAVs entre los instantes de 1.2 segundos y 1.8 segundos. Como se acaba de comentar, mientras que el UAV 1 comienza su trayectoria tras el despegue en el segundo 1.5, el UAV 2 se encuentra realizando una espera hasta que llegue el momento adecuado. 3.1.3 Ejemplo 3 En este tercer ejemplo se va a mostrar una flota de tres UAVs saliendo desde estratopuertos distintos. Al igual que en el ejemplo anterior, esta flota mantendrá una distancia de seguridad entre ellos gracias a las esperas iniciales que realizan. Las tres ubicaciones de los estratopuertos serán aportadas en los datos iniciales, de esta manera se evitan diálogos adicionales. También cabe destacar que se pondrá un tiempo de vuelo bastante superior al resto de tiempos para comprobar que el algoritmo funciona correctamente independientemente del número de vueltas. Los nuevos datos introducidos 3.1 Inspección de un punto 83 son los que se muestran en la Figura 3.12. (a) Código de los datos. (b) Workspace de los datos. Figura 3.12 Datos del ejemplo 3. Como se puede observar en las imágenes, se ha vuelto a mantener la misma velocidad y radio de giro de las aeronaves por simplicidad. A diferencia de los otros ejemplos, en este caso se han añadido dos estratopuertos adicionales ya que se pretende usar una flota de 3 vehículos no tripulados que salgan desde distintas localizaciones. Las posiciones y los ángulos se han elegido para que la trayectoria descrita sea agradable visualmente y que sea fácilmente legible. Los datos correspondientes al campo de visión, a la altitud, y los datos más avanzados como el stepsize o el timestamp se han mantenido con los mismos valores. Finalmente, se ha modificado el tiempo total y el tiempo por vuelta para obtener una trayectoria recorrida por tres UAVs y que realicen la ruta diez veces. A continuación se va a mostrar los cuadros de diálogos resultantes junto con una representación de la ruta creada, Figura 3.13. (a) Mensajes emergentes. (b) Representación. Figura 3.13 Ruta completa del ejemplo 3. En la imagen de la izquierda se observan unos mensajes similares a los ya obtenidos. En estos se vuelve a mencionar cual es el tipo de inspección que se va a utilizar, que en este caso es de punto. Además, a diferencia que en el ejemplo anterior, aunque se vaya a utilizar más de un UAV, no se pregunta si se quieren utilizar los estratopuertos de forma repetida ya que hay el mismo número de estratopuertos que UAVs a utilizar. Si algún proyecto quisiera repetir algún estratopuerto, solo bastaría con introducirlo de forma reiterada en los datos iniciales. Finalmente, se informa por pantalla de que se necesitarán tres HAPS y de que darán un total de diez vueltas para cumplir con las restricciones o condiciones impuestas. En la representación podemos ver las trayectorias correspondientes a los tres UAVs. La primera 84 Capítulo 3. Simulación aeronave despegará desde el punto más al sureste, ya que este punto es el más cercano a la inspección, y el algoritmo elige por proximidad. El segundo UAV saldrá desde el punto más al noroeste. Cabe destacar que aunque la distancia de estos dos estratopuertos al punto de inspección sea la misma, las aeronaves saldrán en tiempos distintos para mantener una distancia de seguridad adecuada. Finalmente, la tercera y última aeronave saldrá desde el estratopuerto situado en el [−1,−1] en el mapa. Por último, en esta ocasión la ruta de inspección o circuito de espera se realiza con un ángulo de 90 grados en sentido antihorario respecto del eje x. A continuación, se van a mostrar una serie de instantes, a lo largo de la simulación, Figura 3.14 y Figura 3.15. Posteriormente se explicará que ocurre en cada imagen. (a) Despegue UAV 1. (b) Aproximación UAVs 1 y 2. (c) Aproximación UAV 3. Figura 3.14 Despegues y aproximación. Ejemplo 3. (a) Inspección equidistanciada. (b) Regreso UAV 1. (c) Regreso UAVs 2 y 3. Figura 3.15 Inspección y aterrizajes. Ejemplo 3. Este ejemplo es muy completo. En la primera figura se pueden distinguir los 3 UAVs situados cada uno en sus puntos de despegue. El UAV 1 está representado en color rojo, el UAV 2 en color verde y el UAV 3 en color azul. Las tres aeronaves se dirigirán hacia el circuito de inspección como siempre, respetando las restricciones de ángulos de entrada y radios de giro gracias a las curvas de Dubins. Primero despegará el UAV rojo y se aproximará al objetivo. Cuando se haya hecho la espera correspondiente, el UAV verde saldrá detrás del rojo para introducirse también en el circuito de inspección. Por último, saldrá el UAV azul de la misma manera que el UAV 2, es decir, respetando la distancia de seguridad. Todo esto se observa en la segunda imagen donde todavía no ha despegado el UAV 3 y los otros dos ya se encuentran en la aproximación. 3.1 Inspección de un punto 85 En la tercera imagen se puede observar como el UAV 1 ya está en el circuito, al igual que el UAV 2, que está alcanzando el objetivo por primera vez. Además, el UAV 3 se encuentra aproximando a una distancia prudencial al objetivo. Posteriormente, en la cuarta gráfica se encuentran los tres UAVs realizando la inspección. Gracias a la distancia que mantienen constante entre ellos, nunca se chocarán. Además, consiguen que siempre haya un UAV viendo el punto de inspección cada 3 segundos, que es el tiempo de inspección introducido. Y, por último, en la quinta y sexta figura la flota de vehículos no tripulados ya ha terminado de dar las 10 vueltas de inspección. A continuación, proceden a volver a sus respectivos estratopuertos para finalizar la ruta completa. Tras explicar como se ha ejecutado la simulación en este ejemplo, se va a realizar una recopilación y análisis de los resultados finales, Figura 3.16 y Figura 3.17, como se ha hecho anteriormente. Figura 3.16 Resultados obtenidos del ejemplo 3 (1). (a) UAV 1. (b) UAV 2. (c) UAV 3. Figura 3.17 Resultados obtenidos del ejemplo 3 (2). En este ejemplo se tiene una estructura con tres filas las cuales corresponden a las trayectorias de cada uno de los UAVs como se puede observar en la Figura 3.16. Esta ruta tiene muchas más coordenadas que los anteriores dos ejemplos. Esto se debe a que en este caso se han dado un total de 10 vueltas al circuito y, por tanto, ha sido necesario almacenar una gran cantidad de puntos. Al igual que en el ejemplo 2, en este caso se tiene que los UAVs 2 y 3 tienen más puntos que el UAV 1. Esto ocurre porque tienen que esperar al momento adecuado para despegar desde el estratopuerto como se observa en la Figura 3.17. En este conjunto de imágenes se muestra el mismo intervalo de tiempo comprendido entre 2.1 y 2.7 segundos. Se observa en las imágenes que, mientras el primer UAV ya está en mitad del vuelo y en movimiento, el segundo UAV está terminando de despegar y comenzando su ruta; y el último UAV todavía no ha a comenzado su despegue. Por este motivo, se tiene una mayor cantidad de puntos en los UAVs que despegan más tarde. 86 Capítulo 3. Simulación 3.1.4 Ejemplo 4 El último de los ejemplos de la inspección del punto servirá para mostrar el potencial del algoritmo. Para esto se añadirán una gran cantidad de aeronaves despegando desde puntos diferentes. El algoritmo tendrá que calcular todas las trayectorias y los tiempos de espera, para que las distancias entre las aeronaves se mantenga siempre igual. En esta ocasión solo se mostrarán los datos utilizados y algunos momentos interesantes de esta ruta. (a) Código de los datos. (b) Workspace de los datos. Figura 3.18 Datos del ejemplo 4. En este caso, Figura 3.18, se ha cambiado por primera vez el radio de giro ya que se va a trabajar en un entorno más grande. La velocidad se ha mantenido la misma y esto se verá reflejado en la simulación. Ya que, las aeronaves al tener que recorrer grandes distancias a una velocidad relativamente pequeña, hace que la simulación tarde más tiempo. Se han situado cinco estratopuertos que se usarán según su proximidad al punto de inspección. Más adelante se observará que se van a necesitar 9 UAVs y, por tanto, estos puntos iniciales no serán suficientes. Además de los estratopuertos, también se a ha modificado el punto objetivo para que haya cierta distancia con las zonas de despegue. Además, el tiempo de inspección se ha disminuido mucho para que puedan introducirse más aeronaves. Esta frecuencia a la que se quiere inspeccionar hace que tengan que estar las aeronaves cercanas para que se pueda cumplir la especificación. A continuación, Figura 3.19, se va a mostrar como quedaría la ruta final junto con los distintos mensajes que se obtienen en la ventana de comandos de MATLAB. (a) Mensajes emergentes. (b) Representación. Figura 3.19 Ruta completa del ejemplo 4. En esta ocasión hay una gran cantidad de trayectorias que se mezclan por causa del gran número 3.1 Inspección de un punto 87 de vehículos en la flota y de los distintos estratopuertos utilizados. Esto se ve reflejado en la imagen de la ruta completa. Además, se puede observar que el número de aeronaves que van a realizar el recorrido es de 9 UAVs como ya se comentó. También se sabía que los puntos de inicio no iban a ser suficientes y, por tanto, se le tiene que indicar al algoritmo que repita los que sean necesarios. Ahora se van a mostrar seis nuevas imágenes, Figura 3.20 y Figura 3.21, donde se observarán distintos momentos de la simulación completa. Se le ha asignado un color diferente a cada UAV para distinguirlos durante la trayectoria. (a) Despegues iniciales. (b) Aproximación UAVs. (c) Comienzo de inspección. Figura 3.20 Despegues y aproximación. Ejemplo 4. (a) Inspección equidistanciada. (b) Regreso UAVs. (c) Aterrizajes finales. Figura 3.21 Inspección y aterrizajes. Ejemplo 4. Se muestra en este conjunto de imágenes la progresión de la simulación a lo largo del tiempo. Todas las aeronaves comienzan en sus respectivos estratopuertos, como se puede observar en la primera imagen. Posteriormente, los vehículos van despegando por zonas manteniendo una distancia de seguridad mientras el resto de la flota espera al instante correcto como se puede ver en la segunda figura. En la tercera figura se muestra como ya han salido todos los vehículos de sus zonas de despegue y algunos ya han comenzado el circuito de inspección. Cabe destacar que las aeronaves no colisionan en ningún momento gracias a las distancias de seguridad que se mantienen entre ellos. En el segundo trío de fotos se observan los siguientes aspectos. En la primera foto se muestra de forma muy representativa la separación que existe entre las aeronaves. Esto permite ser constante con la inspección del punto respetando las condiciones iniciales. Una vez se han dado todas las vueltas correspondientes, que en este caso era una única pasada, la flota de vehículos no tripulados procede a regresar a sus coordenadas iniciales. Esto se puede observar en las dos últimas figuras del conjunto, dónde los UAVs van regresando de manera progresivo en el mismo orden en el que se desplegaron. 88 Capítulo 3. Simulación Finalmente, una vez terminada la simulación se obtienen los resultados de las trayectorias en formato de estructura, Figura 3.22, como en los anteriores ejemplos. Figura 3.22 Resultados obtenidos del ejemplo 4. Como ya se sabía, en esta ocasión ha sido necesario desplegar una flota de 9 vehículos no tripulados. Las trayectorias de cada aeronave queda registrada para su uso en la estructura que se observa en la imagen. Al igual que en los ejemplos anteriores, a medida que el UAV despega más tarde en el tiempo, es necesario que tenga más puntos en la trayectoria que corresponden con las esperas necesarias. Es por ello, que las últimas aeronaves tienen más puntos que las primeras. Se observa una situación curiosa. La trayectoria del UAV 5 es más corta, en número de puntos, que la de las aeronaves 3 y 4. Esto se debe a que aunque el estratopuerto de los UAV 3 y 4 estuvieran más próximos al objetivo, tenían que dar una vuelta más grande debido al ángulo de salida desde el estratopuerto. 3.2 Inspección de un área En esta sección se expondrán varios ejemplos de la misma manera que se hizo con la inspección del punto. Estos ejemplos contendrán los datos utilizados para obtener las trayectorias además de una serie de imágenes que muestren las simulaciones correspondientes. Los ejemplos que se expondrán irán aumentando en dificultad para entender correctamente la finalidad de las trayectorias. En la sección anterior se enseñó el código correspondiente a la simulación de la inspección de un punto. En este segundo tipo de inspección, el algoritmo para conseguir la simulación es el mismo que se encuentra en los Código 3.1 y 3.2 de la sección anterior. De la misma manera se mostrarán las aeronaves representadas como un punto de color en el espacio y el campo de visión como un rectángulo con las dimensiones correspondientes. 3.2.1 Ejemplo 1 Este primer ejemplo de la inspección de área será el más sencillo. En el se mostrará una única aeronave sobrevolando un área rectangular de dimensiones dadas. El UAV únicamente realizará una inspección y volverá al estratopuerto del que despegó. Al igual que en los ejemplos de la inspección de un punto, se omitirá la utilización del ángulo de guiñada debido al propósito de utilizar estas 3.2 Inspección de un área 89 trayectorias en un entorno real. (a) Código de los datos. (b) Workspace de los datos. Figura 3.23 Datos del ejemplo 1. Primera se va a explicar los datos utilizados para este primer ejemplo, Figura 3.23. La velocidad de la aeronave será de 1 metro por segundo a lo largo de todos los ejemplos por comodidad. Al radio de giro de la aeronave se le ha dado un valor de 0.25 metros para que se observe de una manera más clara las trayectorias como se mostrará más adelante. El estratopuerto del que saldrá la aeronave estará situado en las coordenadas [−3,−3.5] y el UAV despegará en dirección Norte, es decir, formando un ángulo de 90 grados con el eje x en sentido antihorario. El campo de visión en este caso será un rectángulo formado por 1 metro de ancho y 0.6 metros de largo. En este método de inspección se hará uso también del factor de solapamiento alpha . Como su propio nombre indica, este parámetro permitirá controlar cuanto porcentaje del campo de visión revisará zonas ya visitadas. La altitud de las aeronaves será constante a lo largo de todas las trayectorias y con valor de 1 metro. La zona de inspección en este caso será un rectángulo con dimensiones de 6×5.4 metros, donde los extremos se encontrarán a las siguientes distancias del origen de coordenadas: el Norte estará a 2 metros del origen, el Sur se encontrará a 2.5 metros del eje x en dirección hacia abajo, el este y el oeste se encuentran ambos a 3 metros del eje y en direcciones contrarias. El tiempo total de vuelo se ha indicado a 40 segundos teniendo en cuenta como se verá en las siguientes imágenes que el tiempo que tarda la aeronave en dar una vuelta es de aproximadamente 39 segundos. En esta inspección de un área, el dato correspondiente a la frecuencia de inspección no se aplica. A continuación se mostrarán otras dos imágenes en la Figura 3.24. La primera corresponderá con los cuadros de diálogos que aparecen en la ventana de comandos de MATLAB cuando se ejecuta el código. La segunda de las imágenes es una representación de la trayectoria sobre la que se va a realizar la simulación. 96 Capítulo 3. Simulación (a) Despegue. (b) Aproximación a la ruta. (c) Inspección (1). Figura 3.36 Despegue, aproximación e inspección. Ejemplo 3. (a) Inspección (2). (b) Inspección (3). (c) Aproximación final. Figura 3.37 Inspección, regreso y aterrizaje. Ejemplo 3. En esta serie de seis gráficas, se muestran distintos momentos de la simulación. En este caso parece que no concuerda el orden pero esto es debido a que se dan dos vueltas en este ejemplo y, por tanto, se vuelven a revisar zonas ya visitadas. La primera imagen corresponde con el despliegue de la flota. En primer lugar despegará la aeronave 1 que recorrerá la ruta azul. Por cercanía, el siguiente vehículo que despegará será el que recorra la ruta amarilla. Por último, el UAV 2 que se encarga de inspeccionar la zona naranja será el último en despegar. En la segunda imagen se observa como los tres UAVs han llegado al inicio de sus respectivas zonas y proceden a inspeccionarlas al mismo tiempo. Esto es muy útil de cara a evitar colisiones en las trayectorias y para mantener a la flota de vehículos organizada. En la tercera imagen se observa como siguen revisando toda la zona sin dejar ningún espacio sin inspeccionar gracias al campo de visión que queda representado encima de los vehículos. En la cuarta imagen se observa como han terminado la primera inspección y las aeronaves se encuentran realizando una última pasada por el inferior de la zona. Aunque en los ejemplos previos se dirigían al final tras esta revisión, en este caso van a dar dos vueltas de inspección. Por tanto, una vez llegado al extremo de la zona objetivo se encauzan a dar otra vuelta como queda representado en la quinta gráfica. Finalmente, tras haber hecho las dos inspecciones como se esperaba, la flota de vehículos no tripulados se dispone a regresar a los puntos de partida y a aterrizar en sus estratopuertos correspondientes. A continuación se van a mostrar los resultados que se obtienen una vez haya terminado la simulación del vuelo, Figura 3.38. 3.2 Inspección de un área 97 Figura 3.38 Resultados obtenidos del ejemplo 3. De la misma manera que en el resto de casos, los resultados quedan guardados en una estructura de tres elementos donde cada uno de ellos representa una de las trayectorias descritas. Cabe destacar que en estos casos el número de coordenadas es muy parecido entre UAVs, esto es debido a que las distancias a los puntos de inicio son muy próximas unas a otras. 3.2.4 Ejemplo 4 Por último, se va a presentar un ejemplo que supera en complejidad a los anteriores. En este último ejemplo se va a ver el potencial que tiene el algoritmo para organizar una flota de vehículos mucho más grande que en los anteriores ejemplos. En este caso se van a incluir unas especificaciones que necesiten de una flota de 12 vehículos para cumplirse. Los datos correspondientes son los que se observan el la Figura 3.39. Figura 3.39 Datos del ejemplo 4. Como se puede observar, se han añadido hasta seis estratopuertos en distintas posiciones. De esta manera se fuerza al algoritmo a que tenga que elegir el estratopuerto más conveniente para cada zona revisada. Además, se ha aumentado el área total a revisar para ver los resultados de una forma más clara y repartida. También se ha modificado el tiempo total de vuelo, este se ha disminuido para que únicamente den una única vuelta de inspección. 98 Capítulo 3. Simulación (a) Mensajes emergentes. (b) Representación. Figura 3.40 Ruta completa del ejemplo 4. Los mensajes que se observan en esta ocasión son muy similares a los de los ejemplos pasados. En estos mensajes se vuelve a mencionar que se va a realizar la inspección de área y que la división de la zona completa se ha hecho en 12 partes diferentes. Como siempre, estas partes serán inspeccionadas cada una por un único UAV. Además, se menciona que no hay estratopuertos suficientes y, por tanto, solicita poder reutilizar alguno de los estratopuertos proporcionados. Por último, se destaca que el tiempo que se tarda en realizar la inspección es de 20.85 segundos. Como se ha introducido un tiempo total de 21 segundos, el algoritmo determina que solo se podrá dar una vuelta al circuito. Por otro lado, en el lado derecho se encuentra una representación del mapa total con sus trayectorias. La forma de las trayectorias es la misma que en los ejemplos anteriores, donde en este caso se han diferenciado con una gran variedad de colores para visualizarlas correctamente. Los estratopuertos vuelven a estar representados como círculos de color azul. Estos se encuentran en los alrededores de la zona objetivo. Se irán eligiendo su uso según proximidad al inicio de cada trozo a inspeccionar. Algunos serán usados por un único UAV y otros puede que hasta por tres aeronaves distintas. Las dos imágenes que se han comentada se encuentran recogidas en la Figura 3.40. A continuación se van a mostrar seis imágenes que corresponden con la simulación de las trayectorias, todas ellas se encuentran en la Figura 3.41 y en la Figura 3.42. (a) Despegue. (b) Aproximación a la ruta. (c) Inspección (1). Figura 3.41 Despegue, aproximación e inspección. Ejemplo 4. 3.2 Inspección de un área 99 (a) Inspección (2). (b) Inspección (3). (c) Aproximación final. Figura 3.42 Inspección, regreso y aterrizaje. Ejemplo 4. La primera de las imágenes muestra la posición de todas las aeronaves antes del despegue. En esta situación habrá algunas que comiencen su movimiento antes y otras después. Esto dependerá de lo cerca o lejos que se sitúe el estratopuerto de la zona a investigar. En la segunda imagen ya se empiezan a acercar algunas de los UAVs a sus zonas de investigación. Se observa que las aeronaves que salen desde estratopuertos más alejados, han tenido que salir con más tiempo ya que la inspección se realizará de forma simultánea para todas. Por eso, se pueden observar algunos vehículos que todavía no han comenzado su vuelo, ya que estarán más cerca de sus respectivas zonas. En la tercera imagen ya se puede observar una situación interesante donde todas las aeronaves se encuentran volando de una forma organizada recorriendo las mismas zonas en los mismos tiempos. Se ha podido comprobar que aunque se incluya una flota con un mayor número de aeronaves, esto no provoca ningún tipo de inconveniente sobre el algoritmo. Sobre la cuarta imagen se sigue observando este vuelo tan sincronizado. La trayectoria está llegando a su fin y las aeronaves se preparan para realizar su último tramo de inspección sobre el borde inferior de la zona a inspeccionar. Esto que se comenta se muestra en la quinta imagen donde todas las aeronaves realizan este tramo que permite regresar a la base mientras se termina de revisar la zona. Finalmente, cada vehículo de la flota regresa a su estratopuerto correspondiente después de haber completado con éxito la misión como se puede ver en la sexta gráfica. Por último, una vez finalizada la simulación, se obtienen los resultados de las trayectorias de todas las aeronaves, Figura 3.43. Figura 3.43 Resultados obtenidos del ejemplo 4. 100 Capítulo 3. Simulación Como en los ejemplos anteriores, estos resultados vienen expresados como una estructura que contiene doce elementos. Cada uno de estos elementos representa una de las aeronaves utilizadas para la inspección de las doce zonas. En estas trayectorias se tienen cada una de las coordenadas en las que tiene que situarse cada aeronave en cada instante de tiempo separadas por un timestamp introducido en los datos. Además de las coordenadas, también se indica la orientación de estas en todo momento. 4 Pruebas Reales En esta última sección del trabajo, se va a realizar una demostración donde una flota de vehículos no tripulados recorra alguna de las trayectorias que se han presentado en ejemplos anteriores o similares. Primero se van a comentar las posibles aplicaciones reales que pueden tener las trayectorias y se hará una breve explicación del entorno donde se van a probar las trayectorias simuladas. Posteriormente se realizarán tres casos distintos donde se mostrará la trayectoria, se realizará la simulación correspondiente y se probarán en un entorno realista. 4.1 Introducción La planificación de trayectorias es un recurso con mucho futuro en el sector agrario que poco a poco se va digitalizando más y más. Además, otros muchos sectores se desarrollarían más haciendo uso de drones que optimicen las inspecciones de sus equipos de trabajo como podría ser una planta de placas solares. Es por ello que a estas trayectorias que se han planteado en este trabajo se les ha buscado una aplicación. También es importante que partiendo de las trayectorias creadas, se pueda realizar una simulación donde se compruebe el comportamiento que tendrá la flota de vehículos que recorrerán estas rutas. Por tanto, se va a utilizar un ejemplo con aplicación real para demostrar la aplicabilidad del algoritmo creado. El ejemplo que se va a tratar se basará en una misión de vigilancia a través de distintos territorios de la comunidad autónoma de Andalucía. Esta misión se dividirá en tres casos de uso distintos. En el primero de los casos se realizará una inspección de un punto situado en los alrededores de la provincia de Córdoba. Para el siguiente de los casos de uso se pretende realizar una misión de inspección sobre todo el territorio de Andalucía utilizando tan solo 3 aeronaves. Por último, en el tercer ejemplo se realizará la misma inspección que se hizo en el segundo pero se podrán utilizar más aeronaves hasta un total de 6 vehículos en la flota. Para probar estas misiones primero habrá que crear las trayectorias que resuelvan el problema de la forma más eficiente. Posteriormente, se hará una simulación para comprobar que las trayectorias cumplen con lo esperado y observar como recorren las trayectorias, como realizan los despegues y como respetan las distancias entre las aeronaves. Finalmente, se terminará el caso probando las trayectorias en una prueba real y comparando los resultados obtenidos con la simulación. Los algoritmos que realizarán estas trayectorias y las simulaciones correspondientes ya han sido explicados a lo largo del trabajo. Es por ello que a continuación se explicará únicamente el entorno real donde se va a trabajar. 101 102 Capítulo 4. Pruebas Reales 4.1.1 Banco de pruebas Para realizar estas pruebas se ha hecho uso de un espacio, Figura 4.1, perteneciente a la empresa CATEC (Centro Avanzado de Tecnologías Aeroespaciales) dónde se pueden realizar pruebas como la cooperación y coordinación de flotas de vehículos [2]. Figura 4.1 Banco de pruebas. Este banco de pruebas dispone de 20 cámaras VICON que forman un sistema de posicionamiento en interiores. Este sistema obtiene la posición y la actitud de los objetos que se encuentre en un volumen de 15x15x5 metros. Este es el motivo por el que las trayectorias creadas en los ejemplos mostrados en capítulos anteriores utilizaban rangos menores a 10 metros de distancia y a unas velocidades bajas, ya que se pretendían realizar pruebas sobre el banco de pruebas mencionado. 4.1.2 Crazyflies Este concepto corresponde con un pequeño dron de cuatro motores que cabe en la palma de la mano, Figura 4.2. Este dron es muy utilizado para realizar pruebas sobre algoritmos en desarrollo, ya que por su pequeño tamaño, su peso tan ligero de tan solo 27 gramos y las velocidad bajas a las que se usan, proporcionan mucha seguridad en los casos en los que alguna parte de los algoritmos en desarrollo fallen. Además, este dron es muy resistente a los golpes y las colisiones a las que pueden ser sometidos en estas pruebas [4]. Figura 4.2 Crazyflies. 4.2 1ª Aplicación 103 Estos drones se utilizarán como una flota para realizar las pruebas que se verán en los siguientes apartados del trabajo. En este trabajo no se va a explicar como se utilizan los mismos ya que se sale de las líneas de trabajo definidas. 4.2 1ª Aplicación 4.2.1 Simulación En esta primera aplicación se va a realizar una inspección sobre una zona que se encuentra al Norte de la ciudad de Córdoba. Esta inspección se realizará con una flota de tres aeronaves que mantengan una distancia de seguridad entre ellos. Esta distancia también permite que se realice la inspección respetando el tiempo que se quiere mantener entre inspección e inspección. Para crear el entorno se ha importado una imagen del mapa de Andalucía. Ha sido necesario modificar la posición de la imagen para que se encuentre adecuadamente orientada a la hora de hacer una figura con ella en MATLAB. Además, se observará en las siguientes imágenes que las trayectorias realizadas están giradas respecto los ejemplos mostrados en la sección anterior. Esto es debido al entorno donde se realizarán las pruebas que tiene su propio sistema de coordenadas y se han adaptado las trayectorias al mismo. Todo esto se irá comentando a lo largo de la sección. Figura 4.3 Trayectoria 1ªaplicación. En la Figura 4.3 se observa el mapa de Andalucía donde se han situado tres estratopuertos distintos bastante distanciados unos de otros. El punto de inspección queda marcado con un cuadrado naranja al Norte de Córdoba en el centro de la imagen. Las rutas que se encuentran representadas en amarillo son las trayectorias que irán recorriendo las aeronaves. Por otro lado, el circuito representado en color verde muestra la trayectoria de inspección que realizara la flota de vehículos no tripulados. A continuación se va a mostrar en la Figura 4.4, como se hizo en las secciones anteriores, una serie de imágenes que corresponderán a momentos distintos a lo largo de la simulación. 104 Capítulo 4. Pruebas Reales (a) Despegue y aproximación. (b) Inspección. (c) Regreso y aterrizaje. Figura 4.4 Simulación 1ªaplicación. En este caso se observan tres imágenes las cuales corresponden a distintos instantes resaltados a lo largo de la simulación. En la primera de las figuras se muestra como la flota se está aproximando en orden al circuito de inspección. La aeronave de que comenzó más al suroeste de la imagen fue la primera en despegar ya que era el UAV que iba a comenzar la inspección y partía desde un puntos bastante lejano. Posteriormente salió el UAV situado más al sur que también parte desde un estratopuerto alejado de la zona de inspección. Por último, desde el estratopuerto situado en la parte superior derecha de la imagen salió el tercer UAV que será el segundo en llegar al punto de inspección por cercanía al mismo. En la imagen central se observa ya a la flota de vehículos no tripulados realizando las inspecciones en el circuito de color verde. En este circuito toda la flota mantiene una distancia de separación igual a un tercio de la distancia total de la ruta, ya que se trata de una flota de tres aeronaves. A lo largo de la simulación se van a realizar 4 vueltas al circuito antes de regresar a los respectivos estratopuertos. Por último, en la tercera imagen la flota de vehículos ya ha realizado todas las vueltas de inspección y procede a regresar a sus puntos de partida y aterrizar en ellos. 4.2.2 Entorno real La última fase de esta aplicación sería probar las trayectorias en una flota de drones reales para verificar que la simulación realizada corresponde con los resultados que se obtendrían en una aplicación real. Para ello se han elegido los drones crazyflies en el banco de pruebas del que se ha hablado en la introducción de esta sección. Gracias a estos dos se han podido probar las trayectorias que se han creado y se ha comprobado que la simulación funciona de forma satisfactoria. A continuación, se van a mostrar tres imágenes, en la Figura 4.5, las cuales corresponden a distintos momentos a lo largo de la prueba real que se grabó en vídeo. (a) Despegue. (b) Inspección. (c) Regreso y aterrizaje. Figura 4.5 Prueba real 1ªaplicación. 4.3 2ª Aplicación 105 Estas tres imágenes muestras momentos cercanos a los representados en la simulación. En la primera imagen la aeronave que se encuentra en la parte superior izquierda es la primera en despegar mientras que las otras están esperando que que llegue su momento para despegar y aproximarse al circuito de inspección. Cabe destacar que como se dijo en los apartados que comentaban las simulaciones, las trayectorias se han creado en base a estas pruebas en las que las aeronaves despegan verticalmente ponderando su posición en cada instante. Una vez terminado el aterrizaje, se puede comprobar que el tiempo transcurrido es el mismo que se indicaba en los resultados de la simulación. En la segunda imagen la flota ya se encuentra realizando la inspección manteniendo la distancia de seguridad correspondiente. En este circuito se realizan 4 vueltas como estaba indicado antes de regresar a los puntos desde los que despegaron. En ningún momento de la prueba los vehículos se salen de la ruta indicada. Finalmente, una vez realizada las vueltas de inspección, proceden a regresar por las rutas específicas para cada UAV a los puntos desde los que despegaron. En la tercera imagen se observa ese momento en el que el primer UAV está aterrizando (a la izquierda de la imagen) y los otros dos están aproximándose a sus estratopuertos correspondientes. Una vez finalizada esta prueba se verifica que los resultados obtenidos en el entorno real corresponden con los que se obtuvieron en la simulación y respetan las trayectorias creadas con el algoritmo correspondiente. 4.3 2ª Aplicación 4.3.1 Simulación En este segundo ejemplo de posibles aplicaciones de las trayectorias se va presentar una serie de rutas para tres UAVs distintos. Estas aeronaves se encargarán de realizar una inspección de toda la superficie de Andalucía. Para conseguir este propósito aprovecharán que disponen de un gran campo de visión que les permite monitorizar gran parte de Andalucía en pocas pasadas. Para disminuir el tiempo total se han introducido un total de 3 vehículos en la flota. Estos vehículos recorrerán Andalucía en dirección de Este a Oeste y viceversa. Además, una vez creadas las trayectorias se realizará una simulación de las mismas con el algoritmo diseñado y por último una demostración real de como se comportarían las aeronaves. Primero se va a mostrar en la Figura 4.6 las distintas rutas que recorrerá cada aeronaves. Estas trayectorias se pueden observar en la siguiente imagen. Figura 4.6 Trayectoria 2ªaplicación. Índice de Figuras 1.1 High Altitude Pseudo Satellites 2 1.2 Histórico de proyectos con HAPS 2 1.3 Tipos de HAPS 3 1.4 Aplicaciones para las trayectorias [10] 4 1.5 Algoritmos planificadores 5 1.6 Algoritmo RRT 5 1.7 Posibles áreas a inspeccionar 6 2.1 Curva de Dubins. Ejemplo básico 9 2.2 Curva de Dubins. Ejemplo de un estacionamiento en coche 10 2.3 Curva de Dubins. Ejemplo de viraje en barco 11 2.4 Curva de Dubins. Ejemplo de giro de 180◦en avión (1) 11 2.5 Curva de Dubins. Ejemplo de giro de 180◦en avión (2) 12 2.6 Posibles curvas de Dubins 12 2.7 Cálculo curva de Dubins RSL (1) 13 2.8 Cálculo curva de Dubins RSL (2) 13 2.9 Cálculo curva de Dubins RSL (3) 14 2.10 Cálculo curva de Dubins RSL (4) 14 2.11 Cálculo curva de Dubins RSL (5) 15 2.12 Cálculo curva de Dubins RSL (6) 16 2.13 Cálculo curva de Dubins RSL (7) 16 2.14 Ejemplo en MATLAB de curva de Dubins 17 2.15 Ejemplo modificado en MATLAB de curvas de Dubins (1) 18 2.16 Ejemplo modificado en MATLAB de curvas de Dubins (2) 19 2.17 Ejemplo del código de Hsin-Yi 20 2.18 Mensaje del código de Hsin-Yi 20 2.19 Campo de visión de una persona 24 2.20 Campo de visión de una aeronave 25 2.21 Curva cerrada con caminos de Dubins 27 2.22 Inspección básica de un punto 29 2.23 Inspección con desvío 29 2.24 Inspección con espera en hipódromo 30 2.25 Ejemplo de patrón de espera 30 2.26 Inspección con espera en forma de infinito 31 2.27 Inspección con espera en hipódromo y tiempo mínimo 32 2.28 Ejemplos puntos auxiliares (1) 36 113 114 Índice de Figuras 2.29 Ejemplos puntos auxiliares (2) 38 2.30 Cálculo de trayectoria respetando timestamp (1) 46 2.31 Cálculo de trayectoria respetando timestamp (2) 46 2.32 Ejemplo de trayectoria genérica 47 2.33 Ejemplo búsqueda en área (1) 48 2.34 Ejemplo búsqueda en área (2) 49 2.35 Otros ejemplos de búsqueda en área 49 2.36 Ejemplo búsqueda en área cerrada 50 2.37 Ejemplo área rectangular dividida 54 2.38 Selección de coordenadas de pasadas de inspección (1) 57 2.39 Selección de coordenadas de pasadas de inspección (2) 58 2.40 Ejemplo de trayectoria genérica (1) 63 2.41 Ejemplo de trayectoria genérica (2) 63 2.42 Aplicación del algoritmo completo (1) 71 2.43 Aplicación del algoritmo completo (2) 71 3.1 Datos del ejemplo 1 77 3.2 Ruta completa del ejemplo 1 78 3.3 Despegue y primera inspección. Ejemplo 1 78 3.4 Segunda inspección y aterrizaje. Ejemplo 1 79 3.5 Resultados obtenidos del ejemplo 1 79 3.6 Datos del ejemplo 2 80 3.7 Ruta completa del ejemplo 2 80 3.8 Despegues y aproximación. Ejemplo 2 81 3.9 Inspección y aterrizajes. Ejemplo 2 81 3.10 Resultados obtenidos del ejemplo 2 (1) 82 3.11 Resultados obtenidos del ejemplo 2 (2) 82 3.12 Datos del ejemplo 3 83 3.13 Ruta completa del ejemplo 3 83 3.14 Despegues y aproximación. Ejemplo 3 84 3.15 Inspección y aterrizajes. Ejemplo 3 84 3.16 Resultados obtenidos del ejemplo 3 (1) 85 3.17 Resultados obtenidos del ejemplo 3 (2) 85 3.18 Datos del ejemplo 4 86 3.19 Ruta completa del ejemplo 4 86 3.20 Despegues y aproximación. Ejemplo 4 87 3.21 Inspección y aterrizajes. Ejemplo 4 87 3.22 Resultados obtenidos del ejemplo 4 88 3.23 Datos del ejemplo 1 89 3.24 Ruta completa del ejemplo 1 90 3.25 Despegue, aproximación e inspección. Ejemplo 1 90 3.26 Inspección, regreso y aterrizaje. Ejemplo 1 91 3.27 Resultados obtenidos del ejemplo 1 91 3.28 Datos del ejemplo 2 92 3.29 Ruta completa del ejemplo 2 92 3.30 Despegue, aproximación e inspección. Ejemplo 2 93 3.31 Inspección, regreso y aterrizaje. Ejemplo 2 93 3.32 Resultados obtenidos del ejemplo 2 (1) 94 3.33 Resultados obtenidos del ejemplo 2 (2) 94 3.34 Datos del ejemplo 3 95 Índice de Figuras 115 3.35 Ruta completa del ejemplo 3 95 3.36 Despegue, aproximación e inspección. Ejemplo 3 96 3.37 Inspección, regreso y aterrizaje. Ejemplo 3 96 3.38 Resultados obtenidos del ejemplo 3 97 3.39 Datos del ejemplo 4 97 3.40 Ruta completa del ejemplo 4 98 3.41 Despegue, aproximación e inspección. Ejemplo 4 98 3.42 Inspección, regreso y aterrizaje. Ejemplo 4 99 3.43 Resultados obtenidos del ejemplo 4 99 4.1 Banco de pruebas 102 4.2 Crazyflies 102 4.3 Trayectoria 1ª aplicación 103 4.4 Simulación 1ª aplicación 104 4.5 Prueba real 1ª aplicación 104 4.6 Trayectoria 2ª aplicación 105 4.7 Simulación 2ª aplicación 106 4.8 Prueba real 2ª aplicación 107 4.9 Trayectoria 3ª aplicación 107 4.10 Simulación 3ª aplicación 108 4.11 Prueba real 3ª aplicación 109 Índice de Códigos 2.1 Función de curva de Dubins de MATLAB 16 2.2 Función de curva de Dubins de MATLAB modificada 18 2.3 Ejemplo del código de Hsin-Yi 20 2.4 Función dubins_curve.m (1) 21 2.5 Función dubins_curve.m (2) 21 2.6 Función dubins_core.m (1) 22 2.7 Función dubins_core.m (2) 22 2.8 Función dubins_core.m (3) 23 2.9 Curva cerrada con caminos de Dubins 27 2.10 Inspección básica de un punto 28 2.11 Cálculo de la trayectoria más corta 35 2.12 Tiempos de la trayectoria más corta 37 2.13 Cálculo del número de UAVs necesarios 37 2.14 Reorganización de zonas de despegue 39 2.15 Cálculo de la trayectoria 2D (1) 39 2.16 Cálculo de la trayectoria 2D (2) 40 2.17 Cálculo de la trayectoria 2D (3) 41 2.18 Cálculo de la trayectoria 2D (4) 41 2.19 Cálculo del momento de despegue para cada UAV 42 2.20 Cálculo de la trayectoria 4D (1) 43 2.21 Cálculo de la trayectoria 4D (2) 44 2.22 Cálculo de la trayectoria 4D (3) 44 2.23 Cálculo de la trayectoria 4D (4) 45 2.24 Representación de la trayectoria 46 2.25 Función rectángulo de inspección 52 2.26 División del área (1) 53 2.27 División del área (2) 54 2.28 Reorganización de zonas de despegue 55 2.29 Cálculo de trayectoria 2D (1) 55 2.30 Cálculo de trayectoria 2D (2) 56 2.31 Cálculo de trayectoria 2D (3) 57 2.32 Cálculo de trayectoria 2D (4) 58 2.33 Cálculo del número de vueltas (1) 58 2.34 Cálculo del número de vueltas (2) 59 2.35 Cálculo de trayectoria 2D (5) 59 117 118 Índice de Códigos 2.36 Cálculo de trayectoria 4D (1) 60 2.37 Cálculo de trayectoria 4D (2) 60 2.38 Cálculo de trayectoria 4D (3) 61 2.39 Cálculo de trayectoria 4D (4) 61 2.40 Representación de la trayectoria 62 2.41 Algoritmo de inspección genérico (1) 64 2.42 Algoritmo de inspección genérico (2) 65 2.43 Algoritmo de inspección genérico (3) 65 2.44 Algoritmo de inspección genérico (4) 66 2.45 Algoritmo de inspección genérico (5) 66 2.46 Algoritmo de inspección genérico (6) 67 2.47 Algoritmo de inspección genérico (7) 67 2.48 Algoritmo de inspección genérico (8) 68 2.49 Algoritmo de inspección genérico (9) 68 2.50 Algoritmo de inspección genérico (10) 69 2.51 Script de datos iniciales (1) 69 2.52 Script de datos iniciales (2) 70 2.53 Código App Designer (1) 72 2.54 Código App Designer (2) 72 3.1 Simulación del UAV 75 3.2 Simulación del campo de visión 76 Bibliografía [1] App designer — The Math Works, Inc.,,https:// es.mathworks.com/ products/ matlab/ appdesigner.html. [2] Banco de pruebas,http:// www.catec.aero/ es/ aviónica-y-sistemas/ equipamiento/ banco-depruebas-o-testbed-en-interiores. [3] Campo visual — wikipedia, la enciclopedia libre,https:// es.wikipedia.org/ w/ index.php?title= Campo_visual&oldid=143021389. [4] Crazyflie 2.1,https:// www.bitcraze.io/ products/ crazyflie-2-1/ . [5] Dubins path — Wikipedia, the free encyclopedia,https:// en.wikipedia.org/ w/ index.php?title= Dubins_path&oldid=1097724697. [6] dubinsconnection — The Math Works, Inc.,,https:// es.mathworks.com/ help/ nav/ ref/ dubinsconnection.html. [7] Estructuras — The Math Works, Inc.,,https:// es.mathworks.com/help/ matlab/ structures.html. [8] Jet stream — Wikipedia, the free encyclopedia,https:// en.wikipedia.org/ w/ index.php?title= Jet_stream&oldid=1106515485. [9] Altran, global leader in innovation ventajas y utilidades futuras de las haps, 2017. [10] European Higher Airspace, Stratobus project by thales alenia space, 2021, https:// higherairspace.eu/ stratobus-project-by-thales-alenia-space. [11] José Ayala, Length minimising bounded curvature paths in homotopy classes, Topology and its Applications 193 (2015), 140–151. [12] José Ayala, David Kirszenblat, and J Hyam Rubinstein, A geometric approach to shortest bounded curvature paths, arXiv preprint arXiv:1403.4899 (2014). [13] Testo SE & Co, Campo de visión, objeto visible más pequeño y zona de medición, (Consultado en 2022). [14] RPAS Drones, Stratobus: una plataforma estratosférica a medio camino entre un dron y un satélite, (2018). [15] Lester E Dubins, On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents, American Journal of mathematics 79 (1957), no. 3, 497–516. [16] Harold H Johnson, An application of the maximum principle to the geometry of plane curves, Proceedings of the American Mathematical Society 44 (1974), no. 2, 432–435. 119 120 Bibliografía [17] Hsin-Yi Kang, A matlab version of dubins curve based on andrew walker’s work, 2016, https:// github.com/ EwingKang/ Dubins-Curve-For-MATLAB. [18] Kevin Chmiela, Dubins path, 2016, https:// www.youtube.com/ watch?v=Nr6EX_0V8wM. [19] Derek Kingston, Steven Rasmussen, and Laura Humphrey, Automated uav tasks for search and surveillance, 2016 IEEE Conference on Control Applications (CCA), IEEE, 2016, pp. 1–8. [20] Steven M LaValle, Planning algorithms, Cambridge university press, 2006. [21] Anna Martí, Haps, los híbridos entre satélite y dron que usará la esa para ampliar la exploración y mejorar las comunicaciones, (2017). [22] International Virtual Aviation Organisation, Perform an holding pattern entry - direct entry, 2022, https:// mediawiki.ivao.aero/ index.php?title=Perform_an_Holding_pattern_entry_-_ direct_entry. [23] José Ricardo Sánchez-Ibáñez, Carlos J Pérez-del Pulgar, and Alfonso García-Cerezo, Path planning for autonomous mobile robots: A review, Sensors 21 (2021), no. 23, 7898. [24] Andrei M Shkel and Vladimir Lumelsky, Classification of the dubins set, Robotics and Autonomous Systems 34 (2001), no. 4, 179–202. [25] ESA Space Solutions, Services enabled by high altitude pseudo satellites (haps) complemented by satellites, (2017). [26] Andrew Walker, Dubins-curves: an open implementation of shortest paths for the forward only car, 2008–, https:// github.com/ AndrewWalker/ Dubins-Curves. [27] Ángel Aller, Aprende a beneficiarte del campo de visión en tus videojuegos, (2021).