Mejora computacional de herramientas de planificación de trayectorias de avión bajo meteorología adversa
Abstract
En este trabajo se plantea la mejora computacional de las herramientas de planificación de trayectorias de avión en presencia de meteorología adversa. El algoritmo que más tiempo consume en las trayectorias de evitación de tormentas es el de generación del grafo de visibilidad. La implementación actual de la que dispone el Grupo de Ingeniería Aeroespacial propone una solución que resulta muy costosa temporalmente cuando el número de nodos del grafo comienza a aumentar. El primer objetivo del proyecto será implementar este algoritmo para, tras ello, programar una rutina avanzada de creación del grafo de visibilidad basada en el método de Lee con el propósito de conseguir una reducción relevante en el tiempo de ejecución del algoritmo con respecto a las implementaciones más sencillas de las que se dispone actualmente.
Full text
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 Mejora computacional de herramientas de planificación de trayectorias de avión bajo meteorología adversa Autor: Narciso Valverde González Tutor: Antonio Franco Espín Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023
Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Mejora computacional de herramientas de planificación de trayectorias de avión bajo meteorología adversa Autor: Narciso Valverde González Tutor: Antonio Franco Espín Profesor Contratado Doctor Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023
Trabajo Fin de Grado: Mejora computacional de herramientas de planificación de trayectorias de avión bajo meteorología adversa Autor: Narciso Valverde González Tutor: Antonio Franco Espín 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 Este trabajo supone la finalización de cuatro años de proceso de crecimiento tanto académico como personal que ha supuesto el estudio del Grado en Ingeniería Aeroespacial. Por ello, me gustaría dedicar mi proyecto a todas aquellas personas sin las que esto no habría sido posible. A toda mi familia, por valorarme mucho más de lo que yo hago. A mis abuelos, por preocuparse siempre por mí, a mi padre y mi hermano por ser apoyos incondicionales y, en especial, a mi madre, quien me guio por el camino de las ciencias y la ingeniería, por ofrecerme su ayuda en todo lo que le es posible y el gran esfuerzo que dedica día tras día para sacarnos a todos adelante. A mis amigos de toda la vida de los grupos de El corral yCaraqueños quienes, a pesar del poco tiempo que les dedico, me regalan unas horas de desconexión y disfrute las cuales me dan la vida y me permiten seguir mis metas. A los amigos que he ido haciendo a lo largo de la carrera, sobre todo, a los pertenecientes al grupo La Zahurda, por el apoyo, la ayuda y los ánimos que nos hemos brindado mutuamente a lo largo de estos cuatro años. A Antonio Franco, por la tutela realizada en este Trabajo Fin de Grado, allanando el camino en todo momento. Su conocimiento, amabilidad y paciencia han permitido generar un clima de trabajo ameno e incrementar mi motivación por el proyecto. A todos los profesores de la escuela que han propiciado que, a lo largo de los cuatro años de carrera, mi pasión acerca de todo lo que rodea al mundo de la aeronáutica se haya incrementado además de mantener mi afán por seguir adquiriendo conocimientos en este ámbito. Narciso Valverde González Grado en Ingeniería Aeroespacial Sevilla, 2023 I
Resumen En este trabajo se plantea la mejora computacional de las herramientas de planificación de trayectorias de avión en presencia de meteorología adversa. El algoritmo que más tiempo consume en las trayectorias de evitación de tormentas es el de generación del grafo de visibilidad. La implementación actual de la que dispone el Grupo de Ingeniería Aeroespacial propone una solución que resulta muy costosa temporalmente cuando el número de nodos del grafo comienza a aumentar. El primer objetivo del proyecto será implementar este algoritmo para, tras ello, programar una rutina avanzada de creación del grafo de visibilidad basada en el método de Lee con el propósito de conseguir una reducción relevante en el tiempo de ejecución del algoritmo con respecto a las implementaciones más sencillas de las que se dispone actualmente. III
XNotación find Comando de MATLAB para encontrar los índices de los elementos no nulos en una matriz o vector unique Comando de MATLAB para obtener valores únicos en las filas o columnas de una matriz o vector rmmissing Comando de MATLAB para eliminar filas o columnas de una matriz que contengan NaN length Comando de MATLAB para obtener la longitud de un vector o la mayor dimensión de una matriz mean Comando de MATLAB para calcular la media de los elementos de una matriz o vector false Comando de MATLAB para crear una matriz cuadrada llena de ceros lógicos polyxpoly Comando de MATLAB para detectar las intersecciones entre dos polilíneas inpolygon Comando de MATLAB para para determinar si un conjunto de puntos se encuentra dentro de un polígono abs Comando de MATLAB para calcular el valor absoluto de un número o los elementos de una matriz graph Comando de MATLAB para la creación de grafos shortestpath Comando de MATLAB para determinar el camino más corto entre dos nodos en un grafo geoplot Comando de MATLAB para trazar gráficos y visualizar datos en coordenadas geográficas en un mapa reducem Comando de MATLAB que reduce el número de puntos de una región según una cierta tolerancia
Índice de Figuras 1.1 Radar meteorológico de un avión de investigación de la NASA (imagen extraída de [1]) 2 1.2 Zona de convergencia intertropical (imagen extraída de [2]) 3 1.3 Restos del vuelo 242 de Southern Airways (imagen extraída de [3]) 4 1.4 Visualización en cabina que ofrece un radar meteorológico (imagen extraída de [4]) 5 2.1 Ejemplos de casos para los que el algoritmo concluirá que los vértices iyjno son visibles 12 2.2 Diagrama de flujo de la primera versión del algoritmo Ingenuo 15 2.3 Recta que une el vértice icon vértice jen una situación genérica 16 2.4 Base no ortogonal formada por los vectores característicos de las aristas adyacentes al vértice i17 2.5 Tipos de vértice en un polígono 17 2.6 Vectores de la base no ortogonal y vector director de la recta que une icon j18 2.7 Descomposición del vectorven función de los vectoresry s18 2.8 Casos en que la recta que une iyjno pasa por el interior de un polígono 19 2.9 Diagrama de flujo de la segunda versión del algoritmo Ingenuo 24 3.1 Trazado del segmento i j 26 3.2 Árbol de búsqueda binario de las aristas intersecadas 26 3.3 Trazado de la semirrecta σ28 3.4 Diagrama de flujo del programa principal del algoritmo de Lee 34 3.5 Casos en que solo existe un vértice por encima del nodo i36 3.6 Diagrama de flujo de la función VISIBLEVERTICESVA(i,V)38 3.7 Caso en que la intersección del segmento σ con los obstáculos poligonales es un vértice y está repetido 41 3.8 Casos en que uno de los puntos de intersección de σcon los obstáculos poligonales en la inicialización de Tes un vértice 42 3.9 Vectoresuy ty distancia al vértice i43 3.10 Caso en que los vértices i,jj0−1yjj0se encuentran alineados 45 3.11 Diagrama de flujo de la inicialización de Ten la función VISIBLEVERTICESVA(i,V)46 3.12 Casos de actualización de T51 3.13 Diagrama de flujo de la comprobación de visibilidad de los vértices y actualización de T en la función VISIBLEVERTICESVA(i,V)54 3.14 Diagrama de flujo de la actualización de Ten la función VISIBLEVERTICESVA(i,V)55 3.15 Intersección de la recta que pasa por iyjcon e157 3.16 Diagrama de flujo de la función VISIBLEVA(i,jj0,e1)59 4.1 Región de meteorología adversa en consideración 61 XI
XII Índice de Figuras 4.2 Envolvente convexa de la región de meteorología adversa en consideración 62 4.3 Diagrama de flujo del algoritmo de Ramer-Douglas-Peucker 63 4.4 Grafo de visibilidad del vuelo Bruselas-Sarajevo 64 4.5 Trayectoria más corta del vuelo Bruselas-Sarajevo 64 4.6 Trayectoria más corta del vuelo Madrid-Kiev 65 4.7 Trayectoria más corta del vuelo Turín-Varsovia 65 4.8 Trayectoria más corta del vuelo Ljubljana-Berna 66 4.9 Trayectoria más corta del vuelo entre O= (43.1◦N,0.730◦W)yD= (44.7◦N,6.75◦E)66 4.10 Trayectoria más corta del vuelo entre O= (44.5◦N,4.20◦E)yD= (50.0◦N,22.5◦E)67 4.11 Representación del ajuste de los puntos obtenidos en la resolución de problemas 68
Índice de Códigos 2.1 Cálculos preliminares de la primera versión del algoritmo Ingenuo 11 2.2 Cómputo del grafo de visibilidad de la primera versión del algoritmo Ingenuo 13 2.3 Obtención del camino más corto en la primera versión del algoritmo Ingenuo 14 2.4 Cálculos preliminares de la segunda versión del algoritmo Ingenuo 19 2.5 Formación de la base no ortogonal centrada en un nodo i en la segunda versión del algoritmo Ingenuo 21 2.6 Comprobación de la visibilidad de los vértices en la segunda versión del algoritmo Ingenuo 22 3.1 Cálculos preliminares del programa principal del algoritmo de Lee 30 3.2 Cálculo de las matrices EyEinv en el programa principal del algoritmo de Lee 32 3.3 Cómputo del grafo de visibilidad en el programa principal del algoritmo de Lee 33 3.4 Operaciones previas a la ordenación de vértices en sentido antihorario de la función VISIBLEVERTICESVA(i,V)34 3.5 Tratamiento del vértice en caso de que solo quede uno tras la aplicación de simetría en la función VISIBLEVERTICESVA(i,V)36 3.6 Ordenación de vértices en sentido antihorario de la función VISIBLEVERTICESVA(i,V)37 3.7 Intersección del segmento σ y establecimiento de conjuntos en la inicialización del conjunto ordenado Ten la función VISIBLEVERTICESVA(i,V)39 3.8 Comprobación de las intersecciones en la inicialización del conjunto ordenado T en la función VISIBLEVERTICESVA(i,V)41 3.9 Almacenamiento de aristas en la inicialización del conjunto ordenado T en la función VISIBLEVERTICESVA(i,V)43 3.10 Ordenación del conjunto ordenado Ten la función VISIBLEVERTICESVA(i,V)45 3.11 Comprobación de visibilidad de los vértices en la función VISIBLEVERTICESVA(i,V)47 3.12 Búsqueda de las aristas almacenadas en Ten la función VISIBLEVERTICESVA(i,V)47 3.13 Caso en que sea necesario eliminar dos aristas en la actualización de T en la función VISIBLEVERTICESVA(i,V)48 3.14 Casos en que sea necesario eliminar una arista e introducir otra en la actualización de T en la función VISIBLEVERTICESVA(i,V)49 3.15 Caso en que es necesario añadir dos aristas en la actualización de T en la función VISIBLEVERTICESVA(i,V)52 3.16 Comprobación de que el segmento i j no pasa por el interior de un polígono del que ies un vértice en la función VISIBLEVA(i,jj0,e1)56 3.17 Comprobación de la intersección de i j con e1en la función VISIBLEVA(i,jj0,e1)58 A.1 Función para el cálculo de distancias loxodrómicas 74 XIII
XIV Índice de Códigos B.1 Primera versión del algoritmo HeapSort 76 B.2 Primera versión de la función auxiliar Heapify del algoritmo HeapSort 77 B.3 Segunda versión del algoritmo HeapSort 79 B.4 Segunda versión de la de la función auxiliar Heapify del algoritmo HeapSort 80 C.1 Algoritmo de Búsqueda Binaria e Inserción 81 D.1 Primera versión del algoritmo Ingenuo 83 D.2 Segunda versión del algoritmo Ingenuo 85 D.3 Programa principal del algoritmo de Lee 87 D.4 Función VISIBLEVERTICESVA(i,V)89 D.5 Función VISIBLEVERTICESOD(i,V)94 D.6 Función VISIBLEVA(i,jj0,e1)98 D.7 Función VISIBLEOD(i,jj0,e1)98
Índice Abreviado Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Códigos XIII Índice Abreviado XV 1 Introducción 1 1.1 Motivación 1 1.2 Objetivo del proyecto 6 1.3 Estructura del documento 7 2 Algoritmo Ingenuo 9 2.1 Introducción al método Ingenuo 9 2.2 Desarrollo e implementación del algoritmo Ingenuo 10 2.3 Mejora del costo computacional del algoritmo Ingenuo 16 3 Algoritmo de Lee 25 3.1 Introducción al método de Lee 25 3.2 Desarrollo e implementación del algoritmo 30 4 Resultados 61 4.1 Región de meteorología adversa en consideración y generación de problemas 61 4.2 Problemas resueltos 63 4.3 Análisis de resultados 67 5 Conclusiones y líneas futuras 71 Apéndice A Cálculo de distancias loxodrómicas 73 Apéndice B Algoritmos de ordenación 75 B.1 Primera versión del método HeapSort 75 B.2 Segunda versión del método HeapSort 79 Apéndice C Algoritmo de Búsqueda Binaria e Inserción 81 XV
XVI Índice Abreviado Apéndice D Códigos de MATLAB 83 D.1 Primera versión del algoritmo Ingenuo 83 D.2 Segunda versión del algoritmo Ingenuo 85 D.3 Algoritmo de Lee 87 Bibliografía 101
Índice Resumen III Abstract V Notación VII Índice de Figuras XI Índice de Códigos XIII Índice Abreviado XV 1 Introducción 1 1.1 Motivación 1 1.2 Objetivo del proyecto 6 1.3 Estructura del documento 7 2 Algoritmo Ingenuo 9 2.1 Introducción al método Ingenuo 9 2.2 Desarrollo e implementación del algoritmo Ingenuo 10 2.2.1 Cálculos preliminares 10 2.2.2 Cómputo del grafo de visibilidad 11 2.2.3 Obtención del camino más corto 13 2.3 Mejora del costo computacional del algoritmo Ingenuo 16 2.3.1 Desarrollo de un criterio basado en las aristas adyacentes al nodo i para determinar si una recta que pasa por este nodo cruza el polígono del que i es vértice 16 2.3.2 Implementación del algoritmo Ingenuo con el criterio desarrollado 19 3 Algoritmo de Lee 25 3.1 Introducción al método de Lee 25 3.2 Desarrollo e implementación del algoritmo 30 3.2.1 Programa principal 30 Cálculos preliminares 30 Cómputo del grafo de visibilidad 33 Obtención del camino más corto 33 3.2.2 Funciones VISIBLEVERTICESVA(i,V)y VISIBLEVERTICESOD(i,V)33 Ordenación de los vértices en sentido antihorario a partir del sentido positivo del eje x33 Inicialización del conjunto ordenado T38 Comprobación de la visibilidad de los vértices y actualización de T45 3.2.3 Funciones VISIBLEVA(i,jj0,e1)y VISIBLEOD(i,jj0,e1)56 XVII
XVIII Índice 4 Resultados 61 4.1 Región de meteorología adversa en consideración y generación de problemas 61 4.2 Problemas resueltos 63 4.2.1 Problema 1. Vuelo Bruselas-Sarajevo 64 4.2.2 Problema 2. Vuelo Madrid-Kiev 65 4.2.3 Problema 3. Vuelo Turín-Varsovia 65 4.2.4 Problema 4. Vuelo Ljubljana-Berna 66 4.2.5 Problema 5. Vuelo entre O= (43.1◦N,0.730◦W)yD= (44.7◦N,6.75◦E)66 4.2.6 Problema 6. Vuelo entre O= (44.5◦N,4.20◦E)yD= (50.0◦N,22.5◦E)67 4.3 Análisis de resultados 67 5 Conclusiones y líneas futuras 71 Apéndice A Cálculo de distancias loxodrómicas 73 Apéndice B Algoritmos de ordenación 75 B.1 Primera versión del método HeapSort 75 B.2 Segunda versión del método HeapSort 79 Apéndice C Algoritmo de Búsqueda Binaria e Inserción 81 Apéndice D Códigos de MATLAB 83 D.1 Primera versión del algoritmo Ingenuo 83 D.2 Segunda versión del algoritmo Ingenuo 85 D.3 Algoritmo de Lee 87 D.3.1 Programa principal 87 D.3.2 Función VISIBLEVERTICESVA(i,V)89 D.3.3 Función VISIBLEVERTICESOD(i,V)94 D.3.4 Función VISIBLEVA(i,jj0,e1)98 D.3.5 Función VISIBLEOD(i,jj0,e1)98 Bibliografía 101
1 Introducción La presencia de regiones de meteorología adversa durante un vuelo compromete su comodidad y el alcance de sus objetivos. En este trabajo, se plantea la implementación de algoritmos basados en la generación del grafo de visibilidad como mecanismos de evitación de estas tormentas. En este capítulo se exponen algunos de los peligros que representan estas zonas y la solución que se propone para evitarlas. Por otra parte, se presentará el objetivo del proyecto, así como la estructura que seguirá este documento. 1.1 Motivación A pesar de que una aeronave no despega si la tripulación no tiene la certeza de poder garantizar la seguridad de pasajeros y mercancías, diversas circunstancias hacen que esta deba enfrentarse con frecuencia a condiciones atmosféricas adversas en la ejecución de sus planes de vuelo. Los pilotos completan un largo y duro entrenamiento para volar en situaciones extremas con difíciles maniobras y, además, adquieren conocimientos acerca de la meteorología, pero la situación climatológica puede sorprenderlos en el transcurso del vuelo y cambiar en un intervalo corto de tiempo. Esto puede provocar problemas en los procesos de aterrizaje y despegue o alteración de los horarios de los vuelos que desembocan en retrasos molestos para los pasajeros. Asimismo, se pueden producir cancelaciones de los trayectos dependiendo del tipo de avión de que se trate y, en el peor de los casos, accidentes de distinta índole. Durante el aterrizaje, una tormenta puede hacer que una aeronave quede atrapada en corrientes y ser movida de un lado a otro por los vientos cruzados, causando que, si este no pudiese llevarse a cabo, el vuelo se desvíe a otro aeropuerto. En el peor de los casos, es posible que la tripulación pierda el control de la aeronave y se produzca un siniestro con consecuencias de diversa importancia. Estos problemas también pueden producirse en el transcurso del despegue, aunque en estos supuestos es posible realizar una cancelación del trayecto. A lo largo del vuelo, una tormenta puede producir muchos efectos incómodos. Los aviones disponen de radares meteorológicos que se sitúan en el morro de la aeronave y pueden recoger señales de las nubes que aparecen delante de la cabina. De esta manera, se intenta realizar el vuelo a una determinada distancia del núcleo de la tormenta (se tiene un margen de seguridad de 20 millas náuticas según [ 8 ]). Sin embargo, estos radares detectan gotas de agua y son ineficaces ante otras 1
2 Algoritmo Ingenuo Una vez entendida la complejidad del problema al que se pretende dar solución, se van a elaborar una serie de algoritmos con tal propósito. En este capítulo, en primer lugar, se va a formular e implementar un algoritmo para la obtención de la trayectoria más corta cuya generación del grafo de visibilidad se base en el método Ingenuo. Seguidamente, se desarrollará y se aplicará un criterio para reducir ligeramente el tiempo de ejecución del mismo sin perder la esencia del método. 2.1 Introducción al método Ingenuo Con el fin de la generación del grafo de visibilidad, para la implementación del método Ingenuo se va a desarrollar un algoritmo exhaustivo que va recorriendo todas las parejas de vértices (i,j) y comprobando si el segmento i j corta a las aristas de los polígonos que conforman el grafo, es decir, a las aristas obstáculo. Este algoritmo viene descrito en el siguiente pseudocódigo ([10]): Pseudocódigo 2.1 VisibilityGraph(V) Entrada: La variable V , formada por una serie de obstáculos poligonales disjuntos definidos por un conjunto ordenado de vértices y dos puntos que no se encuentren en el interior de los polígonos. Salida: Grafo de visibilidad dado en Conex. 1: Inicializar un grafo vacío Conex. 2: para cada vértice i∈V???? hacer 3: para cada vértice j∈[V−i]hacer 4: si el segmento i j no interseca ninguna arista obstáculo ek0entonces 5: Añadir i j aConex. 6: devolver Conex Nótese que V es, en este caso, la concatenación de la variable V presentada en el Pseudocódigo 1.1 y los nodos de origen Oy destino D. Puesto que se designará como n al número de aristas, el número de nodos del grafo será N=n+2 (el número de aristas a los que se añaden los puntos de origen y destino). Si se observa el Pseudocódigo 2.1, la cuarta línea de código representaría una complejidad O(n=N−2) cuya operación dentro del bucle, debido a la tercera línea de código, se repite N−1 veces y, debido a la segunda línea, lo anterior se repite N veces, por lo que el costo computacional total del algoritmo sería de 9
10 Capítulo 2. Algoritmo Ingenuo O[N(N−1)(N−2)] ∼O(N3). En la siguiente sección, se detalla el desarrollo del algoritmo planteado, así como su implementación en MATLAB. 2.2 Desarrollo e implementación del algoritmo Ingenuo En el caso de la implementación del método Ingenuo, se van a combinar los algoritmos descritos en los Pseudocódigos 1.1 y 2.1 en un mismo programa debido a la sencillez del resultado. Como entradas al algoritmo, se tendrán el punto de origen del vuelo, O , el punto de destino, D , y una matriz de vértices de la región tormentosa, V . En esta última estructura se encuentran las coordenadas de los nodos de los distintos polígonos obstáculo del grafo ordenados en sentido horario, de forma que el primer vértice de cada figura se encuentra repetido para cerrar el polígono y, además, los vértices de distintos polígonos se encuentran separados por NaN (Not a Number en MATLAB). Si se denomina Nin al número de puntos dados en V , esta matriz será de dimensiones Nin ×2. La primera coordenada de los puntos O y D será la latitud geográfica y, la segunda, la longitud geográfica. Por su parte, en la primera columna de la matriz V aparecerán las coordenadas de la latitud geográfica de los vértices y en su segunda columna, las coordenadas de la longitud geográfica. Dichos puntos vendrán dadas en grados decimales. 2.2.1 Cálculos preliminares En primer lugar, se introducirán desde el principio los puntos de origen, O , y de destino, D , al final de la matriz V, separados también por NaN. Acto seguido, interesa crear una variable donde únicamente aparezcan los vértices del grafo sin puntos repetidos y sin NaN entre nodos de distintos polígonos. Este parámetro será una matriz que se denominará WPs (Waypoints) y será de dimensiones N×2 . Para hallarla, primeramente, se aplicará el comando unique, pasando también los argumentos de ‘rows’ para indicar que se debe realizar la búsqueda de valores únicos considerando las filas completas de una matriz, y ‘stable’ para indicar que el orden relativo de los valores únicos se debe preservar en el resultado. Una vez calculada esta matriz, se obtendrá el número total de nodos del grafo, N. Tras ello, puesto que se va a crear un grafo plano, es necesario convertir la información meteorológica dada en la Tierra, una región esférica, a un mundo plano. Este objetivo se va a realizar a través de la transformación cartográfica de Mercator, dada por las siguientes ecuaciones ([12]): x=λ−λ0(2.1) y=ln(π 4+φ 2)(2.2) Donde (x,y) son las coordenadas de un punto en el mapa en la proyección de Mercator, λ es la longitud en radianes, φ es la latitud en radianes y λ0 es una longitud media que vendrá dada por la siguiente fórmula:
2.2 Desarrollo e implementación del algoritmo Ingenuo 11 λ0=1 N N ∑ i=1 λi(2.3) Esta transformación se aplicará tanto a las coordenadas dadas en la matriz V , como a las dadas en la matriz WPs, obteniendo como resultado las matrices VMercator yWPs Mercator . En este caso, en la primera columna de las matrices se introducirán las coordenadas en el eje x de los distintos puntos y en la segunda columna, las coordenadas en el eje y. De igual forma, es necesario inicializar la variable que almacene información asociada a la conectividad de los distintos nodos del grafo. Con tal propósito, se empleará una matriz tipo false y será la anteriormente mencionada Conex. Tendrá dimensiones N×N. 1clear all;close all;clc; 2 3%% Insertar puntos de origen y destino (Madrid y Kiev en este caso) 4O = [40.471926, -3.562640]; D = [50.345001, 30.894699]; 5 6%% Insertar polígonos obstáculo 7 8storms = load(‘TFG_Narciso.mat’); 9storms = struct2cell(storms); 10 V = [storms{1} storms{2}]; 11 12 %% CÁLCULOS PRELIMINARES 13 14 V = [V; NaN(1,2); O; NaN(1,2); D]*pi/180; 15 16 % Matriz de WPs y orden del grafo 17 WPs = rmmissing(unique(V, ‘rows’, ‘stable’)); % Primero se suprimen vértices repetidos 18 % con el comando ’unique’ y tras ello se suprimen los NaN con el comando ’rmmissing’ 19 N=length(WPs); 20 21 % Implementación de la transformación de Mercator 22 Lon_0 = mean(WPs(:,2)); 23 x_V = (V(:,2) - Lon_0); y_V = log(tan(pi/4 + V(:,1)/2)); V_Mercator = [x_V, y_V]; 24 x_WPs = (WPs(:,2) - Lon_0); y_WPs = log(tan(pi/4 + WPs(:,1)/2)); WPs_Mercator = ... 25 [x_WPs, y_WPs]; 26 27 % Inicialización de la matriz de adyacencia lógica 28 Conex = false(N,N); 29 % Código 2.1 Cálculos preliminares de la primera versión del algoritmo Ingenuo. 2.2.2 Cómputo del grafo de visibilidad Se recorrerán todos los vértices del grafo y, fijado un vértice, se comprobará la visibilidad con respecto al resto. Se aprovechará el carácter simétrico del problema, por lo que el vértice de origen de la recta, i , viajará a través del intervalo [1,N−1] y, el vértice de destino, j , a lo largo del intervalo [i+1,N]. Fijados los vértices (i,j) , se comprobará la intersección del segmento que los une, i j , con todas las aristas de los obstáculos poligonales. Para ello se hará uso del comando polyxpoly, pasando
12 Capítulo 2. Algoritmo Ingenuo también el argumento ‘unique’ con objeto de evitar intersecciones duplicadas (en el caso de que un punto de intersección fuese un vértice repetido, de esta forma el comando no detectaría dos puntos de intersección exactamente iguales). Tras ello, para decidir acerca de la visibilidad de un vértice, se mirará en primer lugar el número de puntos de intersección. Este, es siempre superior o igual a 2 ya que los propios vértices i y j son detectados como puntos de intersección. Así, si el número de puntos de intersección no es igual a dos, equivaldrá a que la recta que une i y j ha intersecado algún punto que no es ni el vértice i ni el j, por lo que dichos nodos no serán mutuamente visibles. Hay casos en los que, en realidad, aún habiendo más de dos puntos de intersección, los vértices sí serían visibles, como se muestran en la Figura 2.1. No obstante, en estos casos se puede imponer que no existe visibilidad directa entre los vértices i y j , pero sí desde i hasta el primer vértice intersección ( j−1 o w en la Figura 2.1). Así, tras una aplicación recursiva de esta misma idea, aunque j no sea visible desde i , podrá llegarse desde i hasta j recorriendo el mismo camino que si lo fuese. La principal ventaja de este enfoque radica en la disminución en las conexiones del grafo resultante mediante la eliminación de conexiones redundantes. i j−1 j (a) i w j (b) i w w+3 w+4 j (c) Figura 2.1 Ejemplos de casos para los que el algoritmo concluirá que los vértices i y j no son visibles. Lo único que queda por comprobar para conocer si los vértices i y j son visibles es el paso de la recta que los une por el interior de un polígono. Teniendo en cuenta lo ejecutado anteriormente, esto solo puede ocurrir si dichos puntos se encuentran en el mismo polígono. Para ello, se calcula el
2.2 Desarrollo e implementación del algoritmo Ingenuo 13 punto medio del segmento i j y se hace uso del comando inpolygon, que permite determinar tanto si el punto mencionado no se encuentra en el interior de un polígono como si está sobre una de las aristas obstáculo, en cuyos casos i y j serán mutuamente visibles lo cual se indicará en la matriz Conex. 30 %% CÓMPUTO DEL GRAFO DE VISIBILIDAD 31 for i=1:N-1 32 for j=i+1:N 33 % Intersección entre la recta que une i con j y el conjunto de polígonos 34 % obstáculo 35 [xk, yk] = polyxpoly([WPs_Mercator(i,1), WPs_Mercator(j,1)], ... 36 [WPs_Mercator(i,2), WPs_Mercator(j,2)], V_Mercator(:,1), ... 37 V_Mercator(:,2),‘unique’); 38 39 if length(xk) == 2 % Único caso en que se considera que i y j puedan ser 40 % visibles 41 % Cálculo del punto medio 42 midpoint = [(xk(1) + xk(2))/2, (yk(1) + yk(2))/2]; 43 44 % Se comprueba tanto si el punto medio está dentro del polígono como si 45 % está en una arista 46 [in, on] = inpolygon(midpoint(1), midpoint(2), V_Mercator(:,1), ... 47 V_Mercator(:,2)); 48 if in == 0 || (in == 1 && on == 1) 49 Conex(i,j) = 1; Conex(j,i) = 1; 50 end 51 end 52 end 53 end 54 % Código 2.2 Cómputo del grafo de visibilidad de la primera versión del algoritmo Ingenuo. 2.2.3 Obtención del camino más corto Antes del desarrollo de las líneas de código para la obtención del camino más corto que finalizarían la rutina, se recuerda que estas quedan fuera de la complejidad y el tiempo de computación del algoritmo para la generación del grafo de visibilidad. Una vez hallada la matriz Conex , se creará una matriz de adyacencia con pesos que se denominará Adj. Con esa finalidad, se almacenarán tanto los índices (i,j) de los puntos mutuamente visibles (elementos tales que Conexi,j=1 ) como la distancia loxodrómica dlox entre i y j , calculada por medio de una función auxiliar. Esta viene detallada en el Apéndice A. Se debe tener en cuenta que la matriz Adj ha de ser simétrica y, dado que se emplea un conjunto de líneas de código en el que no se aprovecha el carácter simétrico del problema y considerando que el cálculo de las distancias loxodrómicas lleva acarreados ciertos errores de redondeo, se ejecutarán una serie de comandos para reducir el número de cifras significativas característico de las distancias calculadas a 9. Seguidamente, se crea la matriz de adyacencia empleando el comando sparse. Llegado a este punto es necesario hacer un inciso: puesto que se pretende hallar el camino más corto entre O y D , se debería emplear una matriz de adyacencia cuyos pesos fuesen, en lugar de las distancias loxodrómicas, las ortodrómicas, pues estas representan las distancias más cortas entre dos
14 Capítulo 2. Algoritmo Ingenuo puntos en el globo terráqueo. No obstante, puestos que los puntos mutuamente visibles en el grafo se encuentran relativamente cerca, los dos tipos de distancias mencionados serán muy similares entre sí. Posteriormente, se emplea el comando graph para construir el grafo a partir de la matriz Adj y, a continuación, se hace uso del comando shortestpath, el cual implementa el método de Dijkstra , para hallar tanto los vértices a seguir como la distancia recorrida a lo largo del camino más corto. En adición a lo anterior, se proporcionan unos comandos para dibujar tanto el grafo de visibilidad (en este caso, se propone el empleo de un doble bucle for para aprovechar la simetría del problema debido al alto costo temporal del comando geoplot) como la trayectoria más corta de forma independiente. 55 %% OBTENCIÓN DEL CAMINO MÁS CORTO ENTRE O Y D 56 57 % Creación de la matriz de adyacencia con pesos 58 [i_vec, j_vec] = find(Conex); 59 for n_c = 1:length(i_vec) 60 d_lox(n_c) = LoxodromicDistance(WPs(i_vec(n_c),1), WPs(i_vec(n_c),2), ... 61 WPs(j_vec(n_c),1), WPs(j_vec(n_c),2), 0); 62 end 63 64 significant_numbers = 9; % Número de cifras significativas que se considera 65 factor = 10.^(significant_numbers - ceil(log10(abs(d_lox)))); 66 d_lox = round(d_lox .* factor) ./ factor; 67 68 Adj = sparse(i_vec,j_vec,d_lox); 69 70 % Cálculo del camino más corto 71 G = graph(Adj); 72 [path, distance] = shortestpath(G, N - 1, N); 73 74 %% OPERACIONES AUXILIARES AL ALGORITMO 75 76 %%% Para trazar el grafo de visibilidad, descomentar las siguientes líneas de código 77 % figure(1) 78 %fori=1:N 79 % for j = i : N 80 % if Conex(i,j) == 1 81 % geoplot([WPs(i,1), WPs(j,1)]*180/pi, [WPs(i,2), WPs(j,2)]*180/pi); 82 % hold on 83 % end 84 % end 85 % end 86 % 87 %%% Para trazar el camino más corto, descomentar las siguientes líneas de código 88 % figure(2) 89 % geoplot(V(:,1)*180/pi, V(:,2)*180/pi,’b’); 90 % hold on 91 % geoplot(WPs(path,1)*180/pi, WPs(path,2)*180/pi,’r’); 92 % hold on 93 % for i=1:length(path) 94 % geoplot(WPs(path(i),1)*180/pi,WPs(path(i),2)*180/pi,’rd’) 95 % end Código 2.3 Obtención del camino más corto en la primera versión del algoritmo Ingenuo.
2.2 Desarrollo e implementación del algoritmo Ingenuo 15 En la Figura 2.2, se muestra un diagrama de flujo de la rutina desarrollada con el fin de la implementación del algoritmo Ingenuo para la generación del grafo de visibilidad. Este tipo de diagramas introducidos a lo largo del documento presentará la siguiente convención de símbolos: •: representan el inicio y el final del algoritmo. •: representan los datos de entrada que recibe el algoritmo. •: representan procesos u operaciones. •: representan decisiones tras las que son posibles varios caminos alternativos. •: representan líneas de flujo que indican el orden de ejecución de las operaciones. Inicio: primera versión del algoritmo Ingenuo Origen y destino del vuelo Regiones de meteorología adversa Cálculos preliminares y transformación de Mercator ¿Todo par de nodos (i,j) analizado? Intersección de la recta que une el par puntos con las regiones tormentosas ¿El número de puntos de intersección es dos? ¿El punto medio de la recta está en el interior de un polígono? Añadir par de vértices a la lista de visibles Fin: primera versión del algoritmo Ingenuo No No Sí Sí No Sí Figura 2.2 Diagrama de flujo de la primera versión del algoritmo Ingenuo.
16 Capítulo 2. Algoritmo Ingenuo 2.3 Mejora del costo computacional del algoritmo Ingenuo El algoritmo implementado en la sección anterior presenta un gran inconveniente a pesar de su sencillez: el costo temporal del mismo es muy elevado cuando el número de nodos del grafo comienza a aumentar. En esta sección se va a abordar este problema sin perder la esencia del método Ingenuo. Por tanto, se tratará de sustituir el criterio basado en el punto medio para determinar si una arista pasa por el interior de un polígono, ya que conlleva un costo computacional relevante debido a que el comando empleado realiza comprobaciones de un punto con todo el conjunto de nodos del grafo. Para este criterio, se empleará un planteamiento vectorial y, pese a ser usado en este algoritmo cuando el número de puntos de intersección de la recta que une los vértices i y j es dos, se desarrollará de forma que sea válido para cualquier otra situación. 2.3.1 Desarrollo de un criterio basado en las aristas adyacentes al nodo i para determinar si una recta que pasa por este nodo cruza el polígono del que ies vértice Considérese una situación genérica como la mostrada en la Figura 2.3: i i+ i− j Figura 2.3 Recta que une el vértice icon vértice jen una situación genérica. Se va a contemplar una base no ortogonal formada por el vector r , dirigido del vértice i al nodo posterior, i+, y el vector s, orientado del vértice ial nodo anterior, i−:
2.3 Mejora del costo computacional del algoritmo Ingenuo 17 x y r s i i+ i− Figura 2.4 Base no ortogonal formada por los vectores característicos de las aristas adyacentes al vértice i. Para determinar si la recta que une i con j pasa por el interior del polígono del que i es vértice, es necesario, en primer lugar, conocer si este vértice es convexo o cóncavo. Un vértice convexo en un polígono se caracteriza porque el ángulo interior en dicho vértice es inferior a 180◦ (Figura 2.5a). Análogamente, un vértice cóncavo en un polígono se caracteriza porque el ángulo interior en dicho vértice es superior a 180◦(Figura 2.5b). x y r s z i i+ i− (a) Vértice convexo. x y r s z i i+ i− (b) Vértice cóncavo. Figura 2.5 Tipos de vértice en un polígono. Para conocer si un vértice es convexo o cóncavo, sabiendo que los nodos de los distintos polígonos están ordenados en sentido horario, se realiza el producto vectorial de los vectores r y s , cuyo resultado se denominará q . Puesto que estos vectores están contenidos en el plano xy , q tendrá una única componente no nula en el eje z: q=r× s= exeyez rxry0 sxsy0 = 0 0 rxsy−rysx = 0 0 qz (2.4) Con la finalidad de la aplicación del criterio, resulta necesario conocer el orden en que vienen dados los vértices en la matriz V . Como se ha comentado en la sección 2.2, los nodos se dispondrán en sentido horario pero, en caso de que esto no fuese así, se cambiará la organización de los vértices en
24 Capítulo 2. Algoritmo Ingenuo Inicio: segunda versión del algoritmo Ingenuo Origen y destino del vuelo Regiones de meteorología adversa Cálculos preliminares y transformación de Mercator ¿Todo nodo de origen i analizado? ¿Nodo i perteneciente a un polígono? Formación de base no ortogonal y producto vectorial de los vectores de la base ¿Todo nodo de destino j analizado? Intersección de la recta que une i y j con las regiones tormentosas ¿El número de puntos de intersección es dos? ¿Nodo i es el punto de origen? Cálculo de v y aplicación de criterio ¿Se cumple el criterio? Añadir par de vértices a la lista de visibles Fin: segunda versión del algoritmo Ingenuo No Sí Sí No No Sí No Sí No Sí Sí No Figura 2.9 Diagrama de flujo de la segunda versión del algoritmo Ingenuo.
3 Algoritmo de Lee Los códigos de los algoritmos desarrollados en la sección anterior presentan un importante inconveniente: pese a su sencillez, el costo computacional de los mismos cuando el número de nodos de la región comienza a aumentar es muy elevado. Consecuentemente, surge la necesidad de implementar rutinas que reduzcan este tiempo de computación. En este capítulo se va a desarrollar un algoritmo para la generación del grafo de visibilidad basado en el método de Lee, quien expuso la primera solución no trivial al problema. Dicho resultado presenta, teóricamente, un costo computacional significativamente menor que el obtenido por el método Ingenuo. 3.1 Introducción al método de Lee El algoritmo para la generación del grafo de visibilidad va a seguir un enfoque completamente distinto al desarrollado en el capítulo anterior. Es de gran importancia resaltar la influencia que ha tenido la obra de M. de Berg, O. Cheong, M. van Kreveld yM. Overmars [ 13 ] en esta sección. En primer lugar se fijará el vértice i para, tras ello, emplear otras rutinas que comprueben la visibilidad de este mismo vértice con respecto al resto de nodos del grafo. Se seguirá el procedimiento que se describe en el siguiente pseudocódigo: Pseudocódigo 3.1 VisibilityGraph(V) Entrada: La variable V , formada por una serie de obstáculos poligonales disjuntos definidos por un conjunto ordenado de vértices y dos puntos que no se encuentren en el interior de los polígonos. Salida: Grafo de visibilidad dado en Conex. 1: Inicializar un grafo vacío Conex. 2: para todos los vértices i∈Vhacer 3: si ise encuentra en algún polígono obstáculo entonces 4: Conex ←VisibleVerticesVA(i,V) 5: si no 6: Conex ←VisibleVerticesOD(i,V) 7: devolver Conex Nótese que V es, en este caso, la concatenación de la variable V presentada en el Pseudocódigo 1.1 y los nodos de origen Oy destino D. 25
26 Capítulo 3. Algoritmo de Lee En el método Ingenuo se comprobaba la intersección del segmento i j con todas las aristas obstáculo. En el caso del método de Lee, se plantea la posibilidad de usar la información que se obtiene cuando se testea un vértice j , de forma que se pueda comprobar la visibilidad de otros vértices más rápido. Para conseguir esto es necesario, en primer lugar, ordenar los vértices que se encuentran alrededor del nodo i en un sentido cíclico mientras que, a su vez, se va manteniendo información que ayudará a decidir acerca de la visibilidad del siguiente vértice a ser tratado. Considérese una situación genérica como la mostrada en la Figura 3.1. En ella, se traza el segmento que va del nodo i al nodo j , i j . El vértice j será visible si i j no interseca ninguna arista obstáculo antes de llegar a j . Para comprobar esto, en el método de Lee original se propone la construcción de un árbol de búsqueda binario en el que se almacenen las aristas intersecadas por i j denominado T , tal y como se muestra en la Figura 3.2b. Así, es posible comprobar si j se encuentra tras una arista obstáculo, vista desde i. i j Figura 3.1 Trazado del segmento i j. e4 e3 e2 e1 e8 e7 e6 e5 i j (a) Aristas obstáculo intersecadas por el segmento i j. e4 e2 e1 e1e2 e3 e3e4 e6 e5 e5e6 e7 e7e8 (b) Árbol de búsqueda binario de las aristas intersecadas. Figura 3.2 Árbol de búsqueda binario de las aristas intersecadas. A la vez que se tratan los vértices en un orden cíclico alrededor de i se van guardando las aristas obstáculo intersecadas por i j en T. Sus hojas las almacenan en el siguiente orden: 1. En los nodos inferiores, la hoja que se encuentra más a la izquierda almacena la primera arista obstáculo intersecada por i j , la siguiente el segmento que es intersecado tras el primero y así
3.1 Introducción al método de Lee 27 sucesivamente. 2. Los nodos interiores almacenan la arista obstáculo ek0 que se encuentra más a la derecha del subárbol izquierdo, de forma que todas las aristas que se encuentran en el subárbol derecho son intersecadas por i j tras ek0 y todas las que se encuentran en el subárbol izquierdo son intersecadas anteriormente a ek0(excepto la propia ek0). No obstante en este trabajo, en lugar de emplear un árbol de búsqueda binario, se dispondrá el almacenamiento de las aristas en un conjunto ordenado, de forma que los segmentos del árbol mostrado en la Figura 3.2b quedarían en T de acuerdo con 3.1. Así, para encontrar una determinada arista, se empleará un algoritmo de Búsqueda binaria. T= e1 e2 e3 e4 e5 e6 e7 e8 (3.1) Los algoritmos VisibleVerticesVA (i,V) yVisibleVerticesOD (i,V) desarrollan lo descrito anteriormente. La función VisibleVerticesVA (i,V) no tendrá en cuenta los nodos aislados, es decir, se empleará para comprobar la visibilidad de los nodos pertenecientes a polígonos obstáculo alrededor del vértice i cuando este pertenece a algún polígono obstáculo; a través de esta función se halla lo que se denominará "autografo" de visibilidad. La función VisibleVerticesOD(i,V)se empleará para comprobar la visibilidad de todos los nodos alrededor del vértice i cuando este es, o bien el nodo de origen O o bien el de destino D ; se empleará para introducir estos puntos en el "autografo" de visibilidad. Puesto que es necesario realizar una inicialización del conjunto T , se va a trazar una semirrecta denominada σ desde el vértice i en el sentido positivo del eje x como se muestra en la Figura 3.3, de forma que las aristas intersecadas por σse almacenarán en T. Tras lo anterior, se irán visitando los vértices alrededor del nodo i en orden cíclico. Si dos o más vértices se encuentran en la misma dirección vistos desde el nodo i , lógicamente se tomarán primero los vértices cuya distancia a i sea menor pues el planteamiento de la visibilidad se ha de realizar sobre estos vértices de acuerdo con lo expuesto en la sección 2.2.2 y la Figura 2.1. Seguidamente, se ha de determinar la visibilidad de cada vértice. En caso de que no haya arista alguna en T , el vértice jj0 será visible desde i y, en caso contrario, se emplearán las rutinas denominadas VisibleVA (i,jj0,e1) yVisibleOD (i,jj0,e1) que decidirán acerca de la visibilidad del nodo jj0 respecto al vértice i . Finalmente, se actualizará el conjunto ordenado T tratando las aristas incidentes a jj0. El algoritmo VisibleVerticesVA(i,V)viene descrito en el Pseudocódigo 3.2.
28 Capítulo 3. Algoritmo de Lee x y σ i j1 j2 Figura 3.3 Trazado de la semirrecta σ. Pseudocódigo 3.2 VisibleVerticesVA(i,V) Entrada: Una serie de obstáculos poligonales disjuntos definidos por un conjunto ordenado de vértices, V, y un punto ino aislado que no se encuentra en el interior de ningún obstáculo. Salida: El conjunto de todos los nodos visibles desde i. 1: Ordenar los nodos del grafo en sentido cíclico alrededor del vértice i a partir del sentido positivo del eje x . En caso de empates, los nodos más cercanos a i se dispondrán en primer lugar. Sea j1 , j2,...,jn−1la lista ordenada. 2: Sea σ la semirrecta trazada desde i en el sentido positivo del eje x . Encontrar las aristas obstáculo intersecadas por σ y almacenarlas en un conjunto ordenado T en el orden en que son intersecadas por σ. 3: para j0←1an−1hacer 4: si T está vacío entonces 5: Añadir i jj0a Conex 6: si no 7: si VisibleVA(i,jj0,e1)entonces 8: Añadir i jj0aConex. 9: Insertar en Tlas aristas incidentes al vértice jj0que se encuentran en el lado antihorario de la semirrecta σ. 10: Eliminar de Tlas aristas incidentes al vértice jj0que se encuentran en el lado horario de la semirrecta σ. 11: devolver Conex Nótese que Ves, en este caso, la misma variable Vpresentada en el Pseudocódigo 1.1. El algoritmo VisibleVerticesOD (i,V) se rige por el mismo pseudocódigo que el que se muestra para la función VisibleVerticesVA(i,V)con las siguientes salvedades:
3.1 Introducción al método de Lee 29 • En adición a los obstáculos poligonales, recibe como entrada los puntos de origen O y destino D (es decir, en esta función, V será la concatenación de la variable V presentada en el Pseudocódigo 1.1 y los nodos de origen Oy destino D). •El punto isería, en este caso, un nodo aislado. • En lugar de emplear la función VisibleVA (i,jj0,e1) para decidir acerca de la visibilidad de vértices desde el nodo i, se empleará una función denominada VisibleOD(i,jj0,e1). Sumado a ello, las funciones VisibleVerticesVA (i,V) y VisibleVerticesOD(i,V) presentarán otras diferencias en lo referente a su implementación en MATLAB que se detallarán más adelante. La rutina VisibleVA( i,jj0,e1 ) decide si un cierto vértice jj0 es visible desde el nodo i cuando este pertenece a alguno de los polígonos obstáculo. En primer lugar se ha de comprobar si el segmento i jj0 atraviesa el interior del polígono del que i es un vértice, en cuyo caso los nodos no serían mutuamente visibles. En caso contrario, se realiza la intersección de la primera arista almacenada en T , es decir, e1 con el segmento i jj0 . Si no hay intersección o si la intersección es justamente el vértice jj0 , este nodo e i serán mutuamente visibles. Lo explicado se recoge en el Pseudocódigo 3.3. Pseudocódigo 3.3 VisibleVA(i,jj0,e1) Entrada: El vértice jj0sobre el que se plantea la visibilidad. Salida: Verdadero oFalso. 1: si i jj0atraviesa el interior del polígono del que ies un vértice entonces 2: devolver falso 3: si no 4: si i jj0interseca e1y la intersección no es el propio jj0entonces 5: devolver falso 6: si no 7: devolver verdadero Por su parte, la rutina VisibleOD( i,jj0,e1 ) decide si un cierto vértice jj0 es visible desde el nodo i cuando este es el punto de origen O o el de destino D . La diferencia que presenta con respecto a la función VisibleVA( i,jj0,e1 ) radica en que no es necesario comprobar si i jj0 atraviesa el interior del polígono del que ies un vértice, pues en este caso los nodos iestán aislados. En cuanto a la complejidad característica de las funciones VisibleVerticesVA(i,V)y VisibleVerticesOD (i,V) , la primera línea del pseudocódigo 3.2 presentaría una complejidad de O[(n−1)log2(n−1)] ∼O(Nlog2N) para la función VisibleVerticesVA (i,V) y una de O[(N− −1)log2(N−1)] ∼O(Nlog2N) para la función VisibleVerticesOD (i,V) , por lo que la ordenación de vértices presentaría la misma complejidad para ambas funciones. La segunda línea del pseudocódigo 3.2 presenta una de O(n)∼O(N) para las dos funciones, por lo que no se vería incrementada la complejidad total del algoritmo, y tanto la búsqueda en T como su actualización representa, en el peor de los casos, un costo computacional del O(log2n)∼O(log2N)que, puesto que se realiza para todo nodo alrededor del vértice i , lleva a una complejidad de O[(n−1)log2N]∼O(Nlog2N) para la función VisibleVerticesVA (i,V) y una de O[(N−1)log2N]∼O(Nlog2N) para la función VisibleVerticesOD (i,V) . Por ende, la complejidad de la rutinas VisibleVerticesVA (i,V) y VisibleVerticesOD(i,V)es de O(Nlog2N). Puesto que la función VisibleVerticesVA(i,V) se ejecuta para cada vértice perteneciente a polígonos obstáculos, es decir, n veces y VisibleVerticesOD(i,V) se ejecuta para los nodos
30 Capítulo 3. Algoritmo de Lee de inicio y final, es decir, 2 veces, el algoritmo de Lee presentaría una complejidad total de O(nN log2N+2Nlog2N)∼O(N2log2N). En la siguiente sección, se detalla el desarrollo del algoritmo planteado así como su implementación en MATLAB. 3.2 Desarrollo e implementación del algoritmo El algoritmo completo será resultado de la disposición de una serie de códigos: un programa principal que abarcará los algoritmos descritos en los pseudocódigos 1.1 y 3.1 y, por otro lado, las respectivas funciones VisibleVerticesVA (i,V) ,VisibleVerticesOD (i,V) ,VisibleVA (i,jj0,e1) y VisibleOD(i,jj0,e1). A continuación, se desarrollan los algoritmos planteados. 3.2.1 Programa principal Al igual que en el desarrollo del método Ingenuo, como entradas al algoritmo, se tendrán el punto de origen del vuelo, O , el punto de destino, D , y una matriz de vértices de la región tormentosa, V . Cálculos preliminares Se obtiene el vector W , se creará la matriz WPs, se obtendrá el número de nodos N y se implementará la transformación de Mercator de acuerdo a lo descrito en las secciones 2.2.1 y 2.3.2. En adición a lo anterior, se calculará el número de puntos dado en la matriz V , Nin y se computará el número de aristas obstáculo del grafo, n . Además, para la ordenación de los vértices en sentido antihorario descrita en la sección anterior se creará un vector denominado index , el cual almacenará los índices de los puntos del grafo que se ordenarán una vez fijado el nodo i. 1clear all;close all;clc; 2 3%% Insertar puntos de origen y destino (Puntos artificiales en este caso) 4O = [43, -0.86]; D = [44.7, 6.75]; 5 6%% Insertar polígonos obstáculo 7 8storms = load(‘TFG_Narciso.mat’); 9storms = struct2cell(storms); 10 V = [storms{1} storms{2}]; 11 12 %% CÁLCULOS PRELIMINARES 13 14 V = [V; NaN(1,2); O ; NaN(1,2); D]*pi/180; N_in = length(V); 15 16 % Vector que contiene los índices de la matriz V donde se sitúan los NaN’s y obtención 17 % del número de polígonos obstáculo del grafo 18 W=find(isnan(V(:,1))); 19 npol = length(W) - 1; 20 21 % Matriz de WPs, orden del grafo y número de aristas obstáculo del grafo 22 V_0 = V; V_0(W,:) = []; % Primero se suprimen los NaN 23 WPs = unique(V_0, ‘rows’, ‘stable’); % Seguidamente se suprimen los vértices repetidos 24 N=length(WPs); n = N - 2;
3.2 Desarrollo e implementación del algoritmo 31 25 26 % Implementación de la transformación de Mercator 27 Lon_0 = mean(WPs(:,2)); 28 x_V = (V(:,2) - Lon_0); y_V = log(tan(pi/4 + V(:,1)/2)); V_Mercator = [x_V, y_V]; 29 x_WPs = (WPs(:,2) - Lon_0); y_WPs = log(tan(pi/4 + WPs(:,1)/2)); WPs_Mercator = ... 30 [x_WPs, y_WPs]; 31 32 % Vector de índices de los puntos del grafo 33 index = 1:N; 34 % Código 3.1 Cálculos preliminares del programa principal del algoritmo de Lee. Por otra parte, para el desarrollo del algoritmo es interesante crear una variable que, conociendo el índice de una determinada arista obstáculo, informe acerca de cuáles son los vértices extremos de dicha arista. Esta variable se denominará E y será una matriz de dimensiones Nin ×2 . En los casos en los que la pareja de vértices consecutivos incluya un NaN, se considerarán dos NaN en esa fila de E . El empleo de estos índices resulta de interés en la inicialización del conjunto ordenado T cada vez que se llame a la función VisibleVertices (i,V) , ya que, al obtener mediante el comando polyxpoly los índices de los segmentos de polilínea intersecados, estos se corresponden con las filas de la matriz Esegún el planteamiento empleado. Análogamente a lo anterior, es conveniente crear una variable que, conociendo el índice de un determinado vértice, informe acerca de cuáles son las aristas obstáculo incidentes a ese nodo. Esta variable se denominará Einv y será una matriz de dimensiones N×2 . De igual modo, se almacenarán los mismos índices de arista que los tenidos en cuenta para la matriz E. Puesto que la arista k será la contenida entre los vértices j y j+ (vértice posterior a j ), las matrices EyEinv presentarán la siguiente estructura: E= Fila . . . k− k . . . . . .. . . j−j j j+ . . .. . . Einv = Fila . . . j j+ . . . . . .. . . k k− k+k . . .. . . (3.2) Para el cálculo de estas matrices en MATLAB, se hallará el número de polígonos del grafo npol y se irá visitando cada polígono e introduciendo la información pertinente en dichas matrices a partir del vector W . Teniendo en cuenta que p es el índice del polígono del cual se quiere introducir información, el primer vértice de dicho polígono será Wp−1−2(p−1)+1 (a la posición del NaN entre el polígono actual y el anterior se resta el doble del número de polígonos por el que se ha pasado sin contar el actual y se suma una unidad para obtener el primer vértice del polígono) y el último será Wp−2p (a la posición del NaN entre el polígono actual y el siguiente se resta el doble número de polígonos por el que se ha pasado). Por otra parte, la primera arista de un determinado polígono será Wp−1+1 (a la posición del NaN entre el polígono actual y el anterior se suma una unidad) y la última, Wp−2 (a la posición del NaN entre el polígono actual y el siguiente se restan dos unidades). Además, se recuerda que Nvecp es el número de vértices del polígono p.
32 Capítulo 3. Algoritmo de Lee Teniendo en cuenta lo anterior, los elementos de la matriz E presentarán las siguientes expresiones: •Si p=1: Ek,1=ksi k≤Nvec1 Ek,2=k+1si k<Nvec1 1si k=Nvec1 Ek,1=Ek,2=NaN si Nvec1<k≤Nvec1+2 (3.3) •Si 1<p≤npol: Ek,1=k−2(p−1)si Wp−1+1≤k≤Wp−2 Ek,2=k−2(p−1)+1si Wp−1+1≤k<Wp−2 Wp−1−2(p−1)+1si k=Wp−2 Ek,1=Ek,2=NaN si Wp−2<k≤Wp (3.4) Por otra parte, los elementos de la matriz Einv tendrán las siguientes expresiones: •Si p=1: Einvj,1=jsi j≤Wp−2 Einvj,2=Nvec1si j=1 j−1si 1<j≤Wp−2 (3.5) •Si 1<p≤npol: Einvj,1=j+2(p−1)si Wp−1−2(p−1)+1≤j≤Wp−2p Einvj,2=j+2(p−1)+Nvecp−1si j=Wp−1−2(p−1)+1 j+2(p−1)si Wp−1−2(p−1)+1<j≤Wp−2p (3.6) Por último, se inicializará una matriz de adyacencia lógica Conex como se realizó para los métodos Ingenuo. 35 % Cálculo de las matrices E y Einv 36 E=NaN(length(V),2); Einv = zeros(N,2); 37 for p = 1:npol 38 if p == 1 39 E(1: W(p) - 2,:) = [1 : W(p) - 2; 2 : W(p) - 2, 1]’; 40 Einv(1: W(p) - 2,:) = [1 : W(p) - 2; W(p) - 2, 1 : W(p) - 3]’; 41 else 42 E(W(p-1) + 1: W(p) - 2,:) = [W(p-1) + 1 - 2*(p-1) : W(p) - 2*p; ... 43 W(p-1) + 2 - 2*(p-1) : W(p) - 2*p, W(p-1) + 1 - 2*(p-1)]’; 44 Einv(W(p-1) - 2*(p-1) + 1 : W(p) - 2*p,:) = [W(p-1) + 1 : W(p) - 2; ... 45 W(p) - 2, W(p-1) + 1 : W(p) - 3]’; 46 end 47 end 48 49 % Inicialización de la matriz de adyacencia lógica 50 Conex = false(N,N); 51 % Código 3.2 Cálculo de las matrices EyEinv en el programa principal del algoritmo de Lee.
3.2 Desarrollo e implementación del algoritmo 33 Cómputo del grafo de visibilidad Tras las operaciones preliminares, es hora de implementar la generación del grafo de visibilidad. Para ello, se empleará un bucle for en el que se irán fijando los nodos i ; cuando i pertenezca a polígonos obstáculo se llamará a la función VisibleVerticesVA (i,V) , encargada de la generación del "autografo" de visibilidad y, en caso contrario, se llamará a la función VisibleVerticesOD (i,V) , encargada de la introducción de los puntos de origen O y destino D en el "autografo" de visibilidad. 52 %% CÓMPUTO DEL GRAFO DE VISIBILIDAD 53 54 for i = 1:N 55 if i <= n 56 Conex = VisibleVerticesVA(i, V_Mercator(1:N_in - 4,:), WPs_Mercator(1:n,:), ... 57 index(1:n), E, Einv, Conex, N_in); 58 else 59 Conex = VisibleVerticesOD(i, V_Mercator, WPs_Mercator, N, index, E, Einv, ... 60 Conex, N_in); 61 end 62 end 63 % Código 3.3 Cómputo del grafo de visibilidad en el programa principal del algoritmo de Lee. Obtención del camino más corto El proceso seguido para la obtención de la trayectoria más corta es exactamente el mismo que el descrito en la sección 2.2.3 y el Código 2.3. En la Figura 3.4, se muestra un diagrama de flujo del programa principial para la implementación del grafo de visibilidad a través del algoritmo de Lee. 3.2.2 Funciones VISIBLEVERTICESVA(i,V)y VISIBLEVERTICESOD(i,V) En primer lugar, se va a desarrollar la función VisibleVerticesVA (i,V) , empleada para la generación del "autografo" de visibilidad. La función para la introducción de los puntos de origen O y destino Dal "autografo", VisibleVerticesOD(i,V), se programará a partir de la anterior. La implementación de la función se dividirá en tres partes: 1. Ordenación de los vértices en sentido antihorario a partir del sentido positivo del eje x. 2. Inicialización del conjunto ordenado T. 3. Comprobación de la visibilidad de los vértices y actualización de T. Ordenación de los vértices en sentido antihorario a partir del sentido positivo del eje x La primera operación que debe realizar la función es obtener los índices de los puntos que se almacenarán en j (nótese que, en esta ocasión, j no es un único nodo, sino que será un vector de nodos) y se ordenarán alrededor del vértice i.
40 Capítulo 3. Algoritmo de Lee Tras lo anterior, se definirá un índice k0 que recorrerá cada una de las intersecciones por separado, para las que se realizarán una serie de comprobaciones: 1. En primer lugar, se comprueba si el punto de intersección correspondiente a la arista intersecada que se considere no es el propio vértice i , pues las aristas incidentes a este nodo no representan un obstáculo, por lo que no serán almacenadas en T . Se comprueba únicamente la coordenada en el eje x del punto de intersección, pues las coordenadas en el eje y para todos los puntos de intersección son iguales. 2. A continuación, se extrae la arista obstáculo que se ha intersecado de la segunda columna de la matriz K , es decir, la arista k (puesto que el comando polyxpoly detecta las aristas que van del vértice j al j+ , nodo posterior al vértice j , aunque la intersección sea el propio vértice j ) y se comprueba tanto si la intersección mostrada es un vértice como si es un vértice repetido (en caso de dar con un vértice cuyas coordenadas aparecen duplicadas en la matriz V , estas también aparecerán duplicadas en [xk,yk] , en este caso, al recorrer esta última matriz, en la primera pasada se trataría como un vértice y es en la segunda cuando es necesario detectar que está repetido): a) Para comprobar si la intersección es un vértice, puesto que el índice de la arista k se corresponde con el índice de la posición en que se encuentran las coordenadas del nodo j en la matriz VMercator , se comprueba si estas coordenadas son iguales a las del punto de intersección, en cuyo caso se concreta que la intersección es un vértice. b) Como se ha mencionado anteriormente, en el caso de que la intersección sea un vértice repetido, aparecen en los vectores xk e yk las coordenadas de dicho vértice dos veces: en la primera pasada, se detecta que la intersección es un vértice y en la segunda se detecta si dicho vértice está repetido. Para comprobar esto último, en la segunda columna de la matriz K aparecerá una arista de intersección, que, en este único caso, es la que va del vértice repetido j a j− (nodo anterior al vértice j ) y, puesto que en este supuesto ese índice de arista se corresponde con el índice de la posición en que se encuentran las coordenadas del nodo j− en la matriz VMercator , se busca en dicha matriz las coordenadas en la posición de dicha arista pero incrementada en una unidad y se comprueba si estas coordenadas son iguales a las del punto de intersección, en cuyo caso se concreta que la intersección es un vértice repetido (Figura 3.7). En el caso de la intersección de vértices, solo se considerarán las situaciones en que no se da un vértice repetido, pues para ese caso se detectará antes que la intersección es un vértice y se realizarán las operaciones pertinentes.
3.2 Desarrollo e implementación del algoritmo 41 x yσ k−(ken código MATLAB) k ij j− j+ Si jes un vértice repetido, en MATLAB: Primera pasada: VMercator,k= (xkk0,ykk0) Segunda pasada: VMercator,(k+1)= (xkk0,ykk0) Figura 3.7 Caso en que la intersección del segmento σ con los obstáculos poligonales es un vértice y está repetido. 54 % Bucle para tratar las aristas que se deban insertar en T 55 for k0=1:length(xk) 56 57 % Primera comprobación: que el punto de intersección no sea el punto i del que 58 % parte el segmento sigmagorro 59 if xk(k0) ~= WPs_Mercator(i,1) 60 61 % Arista de intersección 62 k = K(k0,2); 63 64 % Segunda comprobación: si la intersección es justamente un vértice o si es 65 % un vértice repetido 66 isvertice = 0; % Variable que se igualará a 1 si se detecta el caso de 67 % vértice 68 isverticerep = 0; % Variable que se igualará a 1 si se detecta el caso de 69 % vértice repetido 70 71 if V_Mercator(k,1) == xk(k0) && V_Mercator(k,2) == yk(k0) 72 isvertice = 1; 73 elseif V_Mercator(k+1,1) == xk(k0) && V_Mercator(k+1,2) == yk(k0) 74 isverticerep = 1; 75 end 76 % Código 3.8 Comprobación de las intersecciones en la inicialización del conjunto ordenado T en la función VisibleVerticesVA(i,V). Siguiendo a las comprobaciones realizadas, se debe decidir qué aristas se han de almacenar en T : • En caso de que la intersección no sea ni un vértice ni un vértice repetido, se almacenará directamente el índice de la arista obstáculo intersecada en T y su distancia al vértice i en ρ . En el conjunto xcomp se introducirá un NaN (pues nunca será necesario mirar este conjunto para ordenar esa arista).
42 Capítulo 3. Algoritmo de Lee •Cuando la intersección sea un vértice, se pueden dar tres casos: – Si las dos aristas obstáculo incidentes al vértice j se encuentran en el lado horario del segmento σ , no se deberá insertar en T ninguna de ellas, pues estas no van a representar un obstáculo (Figura 3.8a). – Si una de las aristas obstáculo incidentes al vértice j se encuentra en el lado horario del segmento σ y la otra se encuentra en el lado antihorario, se deberá insertar en T únicamente esta última pues es la que va a representar un obstáculo (Figura 3.8b). – Si las dos aristas obstáculo incidentes al vértice j se encuentran en el lado antihorario del segmento σ , se deberán insertar ambas en T pues estas van a representar un obstáculo (Figura 3.8c). x y σ k− k ij j− j+ (a) Aristas obstáculo en el lado horario de σ. x y σ k− k ij j− j+ (b) Una arista obstáculo en el lado horario de σ y otra en el lado antihorario. x yσ k− k ij j− j+ (c) Aristas obstáculo en el lado antihorario de σ. Figura 3.8 Casos en que uno de los puntos de intersección de σcon los obstáculos poligonales en la inicialización de Tes un vértice. Para determinar si se debe insertar o no una arista en T , se seguirán los siguientes procedimientos: – En el caso de la arista k− (índice de la arista anterior a k ), se trazará el vector t , dirigido del vértice j al j− , se normalizará (aunque para esta comprobación no sea necesario) siendo el vector resultado t0 y se mirará la componente en el eje y de dicho vector, t0,y . En caso de que esta sea positiva, la arista k−se añadirá a T. Además, se introducirá la distancia del punto de intersección al vértice i en el conjunto ρ y, por si se da el caso de
3.2 Desarrollo e implementación del algoritmo 43 que se deba introducir también la arista k , para desempatar en cuanto a la ordenación de las aristas en T se introducirá la componente en el eje x de t0 , t0,x en el conjunto xcomp . Para proceder en MATLAB, será necesario conocer el índice de la arista obstáculo k− . Puesto que se conoce el índice de la arista k , se computará el de k− del siguiente modo: en primer lugar, se buscará en la matriz E el índice j del nodo extremo de la arista k para luego, buscar en la matriz Einv el índice de la arista k− incidente a dicho nodo j . De esta forma, se tendrán en cuenta los casos degenerados. Tras ello, puesto que el índice k se corresponde con el índice de la posición en que se encuentran las coordenadas del nodo j en la matriz VMercator y el índice k− con el de la posición de las coordenadas del nodo j− , se hallará el vector t a partir de la matriz VMercator en esta ocasión. – En el caso de la arista k , se trazará el vector u , dirigido del vértice j al j+ , y se realizarán operaciones análogas a las ejecutadas en el caso anterior. Nótese que, en este caso, para proceder en MATLAB al cálculo del vector u con la matriz VMercator , no es necesario buscar el índice de la arista k+ pues, aunque este se corresponda con el caso degenerado de que sea la primera arista de un polígono, en la matriz VMercator los primeros vértices de cada polígono vienen repetidos, por lo que se puede computarudirectamente. ∗ Si es necesario introducir tanto k− como k en T , ocupará una posición menor la que presente menor componente en el eje x de su vector característico normalizado. x yσ k− t k u ρk0 ij j− j+ Se introduce k en Tsi uy>0 Se introduce k− en Tsi ty>0 Figura 3.9 Vectoresuy ty distancia al vértice i. 77 if isvertice == 1 78 % Caso en que la intersección es un vértice 79 kminus = Einv(E(k,1),2); 80 81 % Cálculo del vector t 82 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 83 84 % Para incluir la arista o no en T, se observa la componente en el eje 85 % y del vector t
44 Capítulo 3. Algoritmo de Lee 86 if t_0(2) > 0 87 88 % Introducción de la arista correspondiente en T 89 T = [T; kminus]; 90 91 % Cálculo de la distancia de la arista al punto i 92 rho_k0 = xk(k0) - WPs_Mercator(i,1); 93 rho = [rho; rho_k0]; 94 95 % Almacenamiento de la componente x del vector t normalizado por si 96 % se da el caso de que haya que desempatar en cuanto a la 97 % ordenación en T 98 x_comp = [x_comp; t_0(1)]; 99 end 100 101 % Las operaciones con el vector u son análogas a las realizadas con el 102 % vector t 103 u = V_Mercator(k + 1,:) - V_Mercator(k,:); u_0 = u/norm(u); 104 if u_0(2) > 0 105 T = [T; k]; 106 rho_k0 = xk(k0) - WPs_Mercator(i,1); 107 rho = [rho; rho_k0]; 108 x_comp = [x_comp; u_0(1)]; 109 end 110 111 elseif isverticerep ~= 1 112 % La intersección ya no puede ser un vértice 113 T = [T; k]; 114 rho_k0 = xk(k0) - WPs_Mercator(i,1); 115 rho = [rho; rho_k0]; 116 x_comp = [x_comp; NaN]; % Puesto que en este caso no será necesario 117 % desempatar en cuanto a la ordenación en T, se introduce un NaN 118 end 119 end 120 end 121 % Código 3.9 Almacenamiento de aristas en la inicialización del conjunto ordenado T en la función VisibleVerticesVA(i,V). Posteriormente a la inserción de las aristas en T , es necesario realizar la ordenación de las mismas, donde irán primero las que presenten menor distancia al nodo i y, en caso de empates, se mirará en el conjunto xcomp las componentes en el eje x de los vectores característicos de las aristas normalizados. Para la ordenación, se creará una matriz TreeData donde se introduzcan los parámetros nombrados con el fin de volver a emplear el método HeapSort pero, en este caso, puesto que las condiciones de permutación son distintas a las consideradas para la ordenación de vértices en sentido antihorario a partir del sentido positivo del eje x , se empleará una segunda versión de la función donde se implementen los nuevos supuestos de permutación. Esta se muestra en el Apéndice B.2. TreeData = e1ρ1xcomp,1 e2ρ2xcomp,2 . . .. . .. . . enTρnTxcomp,nT (3.9)
3.2 Desarrollo e implementación del algoritmo 45 Téngase en cuenta que, si el número total de aristas en T es nT , la ordenación en este caso supone un costo computacional de O(nTlog2nT) no considerado anteriormente pero, ni en el peor de los escenarios, esto supone un aumento de la complejidad del algoritmo completo. Por último, tras realizar la ordenación, se está en disposición de rellenar el conjunto Tinv con la posición de las aristas que aparecen en T, finalizando esta etapa del algoritmo. 122 TreeData = [T, x_comp, rho]; 123 124 % Ordenación a través del método HeapSort 125 TreeData = HeapSort_v2(TreeData); 126 T = TreeData(:,1); rho = TreeData(:,3); 127 128 % Relleno del conjunto Tinv 129 for y = 1:length(T) 130 Tinv(T(y)) = y; 131 end 132 % Código 3.10 Ordenación del conjunto ordenado Ten la función VisibleVerticesVA(i,V). En la Figura 3.11, se muestra un diagrama de flujo de las operaciones a realizar para la inicialización del conjunto T en la función VisibleVerticesVA(i,V). Comprobación de la visibilidad de los vértices y actualización de T Se comprobará la visibilidad de los vértices en el supuesto de que j0=1 o en caso de que j0>1 y el segmento i jj0 no se encuentre en la misma dirección que el segmento i jj0−1 , de forma que no se considerarán las situaciones redundantes expuestas en la sección 2.2.2 y la Figura 2.1. En este supuesto, se tiene la situación genérica presentada en la Figura 3.10: x y i jj0−1 jj0 Figura 3.10 Caso en que los vértices i,jj0−1yjj0se encuentran alineados. Para obtener una condición que implique que dichos vértices están alineados, se va a considerar el vector dirigido de i a jj0−1 , # » i jj0−1 y el vector director del segmento i j , v . Si estos vectores van en la misma dirección, el producto vectorial # » i jj0−1×vserá nulo: # » i jj0−1×v= exeyez i jj0−1,xi jj0−1,y0 vxvy0 = 0 0 i jj0−1,xvy−i jj0−1,yvx = 0 0 0 (3.10)
46 Capítulo 3. Algoritmo de Lee Inicio: inicialización de T Intersección del segmento σ con todos los polígonos obstáculo ¿Todos los puntos de intersección analizados? ¿El punto de intersección es el propio vértice i ? ¿El punto de intersección es un vértice? Comprobaciones e introducción en T de las aristas pertinentes ¿El punto de intersección es un vértice repetido? Introducción en T de la arista en que se encuentra el punto de intersección Ordenación de las aristas introducidas en T Fin: inicialización de T No Sí No Sí Sí No No Sí Figura 3.11 Diagrama de flujo de la inicialización de Ten la función VisibleVerticesVA(i,V). Por tanto, se obtiene que la condición para que los segmentos i jj0 y i jj0−1 no se encuentren alineados será pues i jj0−1,xvy−i jj0−1,yvx=0.
3.2 Desarrollo e implementación del algoritmo 47 Tras lo anterior, se realizará la comprobación de si T está vacío, en cuyo caso i y jj0 serán mutuamente visibles, lo cual se indicará a través de la variable bandera vis nuevamente. En caso contrario, la visibilidad de dichos vértices se realizará a través de la función VisibleVA (i,jj0,e1) , que recibirá los argumentos necesarios para su funcionamiento. En caso de obtener vértices mutuamente visibles, se almacenará dicha información en la matriz Conex. 133 %% COMPROBACIÓN DE LA VISIBILIDAD DE LOS VÉRTICES 134 135 for j0 = 1:length(j) 136 if j0 == 1 || (j0 > 1 && (WPs_Mercator(j(j0-1),1) - WPs_Mercator(i,1))* ... 137 (WPs_Mercator(j(j0),2) - WPs_Mercator(i,2)) ~= ... 138 (WPs_Mercator(j(j0),1) - WPs_Mercator(i,1))* ... 139 (WPs_Mercator(j(j0-1),2) - WPs_Mercator(i,2))) % No se considerarán 140 % casos redundantes 141 142 if isempty(T) 143 vis = 1; 144 else 145 vis = VisibleVA(j(j0), i, WPs_Mercator, T(1), E, Einv); 146 end 147 if vis == 1 148 Conex(i,j(j0)) = 1; Conex(j(j0),i) = 1; 149 end 150 end 151 % Código 3.11 Comprobación de visibilidad de los vértices en la función VisibleVerticesVA (i,V) . Por último, se realizará la actualización de T . Se comprobará que el segmento i jj0 no vaya dirigido según el sentido positivo del eje x , pues en caso contrario no hay nada que actualizar porque es lo que se ha realizado en la inicialización de T a través de la intersección del segmento σ con todos los polígonos. La condición que no se ha de cumplir para que i jj0 no vaya dirigido según el sentido positivo del eje xes jj0,y=iyyjj0,x>ix. Por otra parte, también se testeará si el vértice jj0 no es el último en la lista ordenada ya que, en este supuesto, no es necesario actualizar T debido a que no se va a comprobar la visibilidad de más vértices. Acto seguido, se buscarán en la matriz Einv cuáles son las aristas incidentes al vértice jj0 , k y k− , así como la posición que ocupan dichas aristas en el árbol T , posk y posk− (en caso de que las aristas no aparezcan en T, estos valores serán NaN). 152 if ~(WPs_Mercator(j(j0), 2) == WPs_Mercator(i,2) && WPs_Mercator(j(j0),1) ... 153 > WPs_Mercator(i,1)) && j0 < length(j) % Para alpha=0 ya se ha 154 % realizado la inicialización de T y para el último vértice no es 155 % necesario actualizar % T pues no se comprobará la visibilidad 156 % de ningún otro nodo 157 158 %% ACTUALIZACIÓN DE T 159 160 % Búsqueda de la existencia en T tanto de la arista k como de la arista 161 % k^{-} 162 k = Einv(j(j0),1); kminus = Einv(j(j0),2);
48 Capítulo 3. Algoritmo de Lee 163 pos_k = Tinv(k); pos_kminus = Tinv(kminus); 164 % Código 3.12 Búsqueda de las aristas almacenadas en Ten la función VisibleVerticesVA(i,V). A partir de ahora, se explicarán los distintos casos posibles de actualización de T , ordenados de menor a mayor dificultad: •Caso en que sea necesario eliminar ambas aristas en T(Figura 3.12a) En este supuesto, posk= NaN y posk−= NaN. En primer lugar, se cambiará la posición que ocupan en el conjunto T , la cual viene dada en Tinv, por un NaN. Tras ello, será necesario realizar una corrección en Tinv de la posición de las aristas situadas tras k y k− en T , a la que habrá que restar dos unidades (se suprimen dos aristas) y, por último, se eliminarán las aristas de T así como sus distancias al punto i de ρ . En 3.11 se muestra un ejemplo de estas operaciones: Tinv = Fila . . . k k− . . . ek0 . . . . . . posk=NaN posk−=NaN . . . posek0=posek0−2 . . . T= Fila . . . posk− posk . . . posek0=posek0−2 . . . . . . k− k . . . ek0 . . . ;(3.11) 165 % Caso en que sea necesario eliminar ambas aristas en T 166 if ~isnan(pos_k) && ~isnan(pos_kminus) 167 168 Tinv([k; kminus]) = NaN; 169 if pos_kminus > pos_k 170 Tinv(T(pos_kminus + 1:length(T))) = Tinv(T(pos_kminus + ... 171 1:length(T))) - 2; 172 else 173 Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1:length(T))) - 2; 174 end 175 T([pos_k; pos_kminus]) = []; 176 rho([pos_k; pos_kminus]) = []; 177 % Código 3.13 Caso en que sea necesario eliminar dos aristas en la actualización de T en la función VisibleVerticesVA(i,V). •Caso en que sea necesario eliminar una arista e introducir otra (Figura 3.12b) – Si posk= NaN y posk−= NaN, se deberá eliminar de T la arista k y añadir la arista k− en la misma posición en que se encuentra la arista k , eso sí, siempre que k− no se encuentre en la misma dirección que la línea de escaneo, pues en este caso no representaría un obstáculo. Con el objeto de comprobar esto, se calcularán el vector director normalizado del segmento i jj0 , v0 , y el vector normalizado dirigido del vértice j al j− , t0 , y se comprobará si v0= t0o bien si v0=− t0.
3.2 Desarrollo e implementación del algoritmo 49 ∗ Si se cumple lo anterior, se cambiará la posición que ocupa la arista k en el conjunto T dada en Tinv por un NaN. Tras ello, será necesario realizar una corrección en Tinv de la posición de las aristas situadas tras k en T , a la que se requerirá restarle una unidad (se suprime una arista) y, por último, se eliminará la arista k de T y su distancia al punto ide ρ. En 3.12 se muestra un ejemplo de estas operaciones: Tinv = Fila . . . k . . . ek0 . . . . . . posk=NaN . . . posek0=posek0−1 . . . T= Fila . . . posk . . . posek0=posek0−1 . . . . . . k . . . ek0 . . . ;(3.12) ∗ Si lo anterior no se da, se cambiará la posición que ocupa la arista k en el conjunto T dada en Tinv por un NaN y se modificará la posición que ocupa la arista k− en el conjunto T dada en Tinv por la que ocupaba la arista k , posk . Por último, se cambiará en T la arista k por la arista k− . En 3.13 se muestra un ejemplo de estas operaciones: Tinv = Fila . . . k k− . . . ek0 . . . . . . posk=NaN NaN =posk . . . posek0 . . . T= Fila . . . posk . . . posek0−1 . . . . . . k=k− . . . ek0 . . . ;(3.13) –Si posk=NaN yposk−=NaN se daría un caso análogo al anterior. 178 % Casos en que sea necesario eliminar una arista e introducir otra (si esta 179 % no va en la dirección del segmento ij_{j0}) 180 elseif ~isnan(pos_k) && isnan(pos_kminus) 181 182 % Vector director del segmento ij_{j0} 183 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 184 185 % Vector t 186 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 187 188 if (abs(t_0(1) - v_0(1)) < 1e-12 && abs(t_0(2) - v_0(2)) < 1e-12) ... 189 || (abs(t_0(1) + v_0(1)) < 1e-12 && abs(t_0(2) + v_0(2)) ... 190 < 1e-12) % Se ha aplicado una cierta tolerancia 191 % Caso en que solo se suprime una arista 192 Tinv(k) = NaN; Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1:... 193 length(T))) - 1; 194 T(pos_k) = []; 195 rho(pos_k) = []; 196 else
56 Capítulo 3. Algoritmo de Lee 3.2.3 Funciones VISIBLEVA(i,jj0,e1)y VISIBLEOD(i,jj0,e1) Para cerrar el algoritmo de Lee, se va a desarrollar la función VisibleVA (i,jj0,e1) , empleada para comprobar la visibilidad de puntos pertenecientes a polígonos obstáculo. La rutina para testear la visibilidad de los nodos del grafo desde los puntos de origen O y destino D se programará a partir de la anterior. Debido a las entradas que recibe esta función, j denotará, en este caso, un único vértice y no un conjunto como en las funciones VisibleVerticesVA(i,V)yVisibleVerticesOD(i,V). Esta rutina, en primer lugar, creará una variable bandera denominada vis , que será igual a 1 (verdadero) si los vértices son mutuamente visibles y 0 (falso) si no lo son. Seguidamente, se deberá comprobar que el segmento i j no pasa por el interior de un polígono del que i es un vértice. Para ello, se volverá a aplicar el criterio basado en las aristas adyacentes al vértice i que se desarrolló en la sección 2.3.1. Los vectores r y s se computarán del mismo modo que se realizó para la función VisibleVerticesVA(i,V). 1function vis = VisibleVA(j, i, WPs_Mercator, e_1, E, Einv) 2 3vis = 0; 4 5% Comprobación de que el segmento ij no pase por el interior de un polígono del que 6% i es vértice 7 8% Cálculo del vector v 9v = (WPs_Mercator(j,:) - WPs_Mercator(i,:))’; 10 11 % Formación de la base no ortogonal de vectores 12 r = (WPs_Mercator(E(Einv(i, 1), 2),:) - WPs_Mercator(i,:))’; 13 s = (WPs_Mercator(E(Einv(i, 2), 1),:) - WPs_Mercator(i,:))’; 14 15 % Producto vectorial de r y s 16 q_z = r(1)*s(2) - r(2)*s(1); 17 18 % Resolución del sistema de ecuaciones 19 a = (v(1)*s(2) - v(2)*s(1))/q_z; 20 b = (r(1)*v(2) - r(2)*v(1))/q_z; 21 22 % Aplicación del criterio con una cierta tolerancia 23 if (sign(q_z) == -1 && ~(a>1e-12 && b>1e-12)) || (sign(q_z) == 1 && (a>=-1e-12 && ... 24 b>=-1e-12)) 25 % Código 3.16 Comprobación de que el segmento i j no pasa por el interior de un polígono del que i es un vértice en la función VisibleVA(i,jj0,e1). Si lo anterior se cumple, se comprobará si existe la intersección del segmento i j con la primera arista del conjunto T , e1 . Para ello, se empleará el siguiente procedimiento: considérese la situación genérica mostrada en la Figura 3.15:
3.2 Desarrollo e implementación del algoritmo 57 x y AB iA AP v iP e1 A P B i j Figura 3.15 Intersección de la recta que pasa por iyjcon e1. Sea P el punto de intersección de la recta que pasa por los vértices i y j con e1 y sean A y B los nodos extremos de e1 . Sea el iP el vector dirigido del nodo i a P , iA el dirigido i al punto A y AP del dirigido de Aal nodo P. El vector iP se puede expresar del siguiente modo: iP = iA + AP (3.15) De 3.15 se conoce iA (pues se conocen las coordenadas de los nodos i y A ), pero tanto iP como AP son desconocidos. No obstante, son dato tanto v (recuérdese que es el vector dirigido de i a j ) como AB, vector dirigido de AaB, por lo que es posible expresar: iP =βv; AP =γ AB (3.16) Donde β y γ son dos escalares que, a priori, no son conocidos. Sin embargo, introduciendo 3.16 en 3.15 y desarrollando: βv= iA +γ AB −→ βv−γ AB = iA −→ vx−ABx vy−AByβ γ=iAx iAy(3.17) Se obtiene un sistema del que se hallará βempleando la regla de Cramer: β= iAx−ABx iAy−ABy vx−ABx vy−ABy =−iAxABy+iAyABx −vxABy+vyABx (3.18) Con el objetivo de calcular los vectores iA , e iB , se buscarán los índices de los puntos A y B entrando con el índice de la arista obstáculo e1 en la matriz E , para luego buscar las coordenadas de los nodos en la matriz WPsMercator. A fin de comprobar entonces tanto si hay intersección como si la hay pero es justamente uno de los vértices A o B , la norma del vector v , ∥v∥ , ha de ser menor o igual que la del vector iP , iP , lo cual lleva a lo siguiente: ∥v∥≤ iP =∥βv∥=β∥v∥−→ 1≤β(3.19)
58 Capítulo 3. Algoritmo de Lee Si se cumple la condición anterior, iyjserán mutuamente visibles. 26 % Vectores iA y AB 27 iA = WPs_Mercator(E(e_1,1),:)’ - WPs_Mercator(i,:)’; 28 AB = WPs_Mercator(E(e_1,2),:)’ - WPs_Mercator(E(e_1,1),:)’; 29 30 % Cálculo de beta 31 beta = (-iA(1)*AB(2) + iA(2)*AB(1))/(-v(1)*AB(2) + v(2)*AB(1)); 32 33 % Aplicación del criterio con una cierta tolerancia 34 if 1 <= beta + 1e-10 35 vis = 1; 36 end 37 end 38 end Código 3.17 Comprobación de la intersección de i j con e1en la función VisibleVA(i,jj0,e1). En la Figura 3.16, se muestra un diagrama de flujo de la función VisibleVA(i,jj0,e1). Por otra parte, el código desarrollado para la función VisibleOD (i,jj0,e1) resulta muy similar al desarrollado para la función VisibleVA (i,jj0,e1) con una diferencia fundamental: no es necesario comprobar si el segmento i j pasa por el interior de un polígono del que i es vértice, pues los puntos de origen Oy destino Dno están contenidos en ningún polígono. Puesto que la función VisibleOD (i,jj0,e1) es muy semejante a la que se acaba de desarrollar, el código MATLAB implementado se adjunta en el Apéndice D.3.5.
3.2 Desarrollo e implementación del algoritmo 59 Inicio: función VisibleVA(i,jj0,e1) Nodos i y j , e1 y matrices necesarias para realizar cálculos Aplicación del criterio basado en las aristas adyacentes al vértice i ¿ i j atraviesa un polígono del que i es vértice Planteamiento vectorial y cálculo de β ¿ β≥1 ? i y j mutuamente visibles (Verdadero) i y j no son mutuamente visibles (Falso) Fin: función VisibleVA(i,jj0,e1) No Sí Sí No Figura 3.16 Diagrama de flujo de la función VisibleVA(i,jj0,e1).
4 Resultados Tras el desarrollo de códigos basados tanto en el método Ingenuo como en el método de Lee, es necesario comprobar su correcto funcionamiento y obtener una serie de resultados con el fin de analizar la mejoría temporal del algoritmo de Lee con respecto al Ingenuo. En este capítulo se mostrará, en primer lugar, la región de meteorología adversa que se considerará en este trabajo. Tras ello, se seleccionarán puntos de origen y destino para la generación de problemas y se reducirá el número de puntos de la región tormentosa con el propósito de analizar la dependencia del tiempo de ejecución de los algoritmos con el mencionado número de puntos. 4.1 Región de meteorología adversa en consideración y generación de problemas La región tormentosa que se considerará en este trabajo se muestra en la Figura 4.1: Figura 4.1 Región de meteorología adversa en consideración. Si se traza la envolvente convexa de la región, es observable que las tormentas constituyen una franja de polígonos en Europa. 61
62 Capítulo 4. Resultados Figura 4.2 Envolvente convexa de la región de meteorología adversa en consideración. Dicha región está compuesta por un total de n=1196 puntos, a los que será necesario añadir los puntos de origen O , y destino D . Para la selección de estos puntos, se considerarán vuelos entre grandes ciudades de Europa que se encuentren tanto en la proximidad de la envolvente convexa de la región como en el interior cuyas coordenadas geográficas se han tomado de [ 17 ]. Por otra parte, también se ejecutarán trayectos entre puntos artificiales para los que la determinación de la trayectoria más corta presente una dificultad superior. Cada uno de los problemas considerados se resolverá tanto a través de los códigos desarrollados según el método Ingenuo como mediante el código implementado siguiendo el método de Lee, comprobando siempre que los resultados obtenidos para cada problema, empleando los distintos códigos, son idénticos. Por otro lado, será necesario analizar los tiempos de ejecución empleados por cada programa. Además, resulta interesante obtener la sensibilidad de los tiempos de ejecución con respecto al número de nodos del grafo en consideración para analizar empíricamente la complejidad de los métodos y comprobar si se corresponde con los modelos teóricos. Con este fin, se reducirá el número de puntos de la región de meteorología adversa mostrada en la Figura 4.1 empleando el algoritmo de Ramer-Douglas-Peucker ([ 18 ]). Este algoritmo trata de reducir el número de vértices existentes en una polilínea. Se basa en el cálculo de puntos críticos a partir de una tolerancia lineal los cuales constituirán la línea simplificada. Estos puntos, serán los que vayan alcanzando una distancia perpendicular superior a la tolerancia impuesta, respecto a la línea base a tener en cuenta que, en principio, es la limitada por los puntos inicial y final. Después, se calcularán las distancias perpendiculares de todos los puntos intermedios: • Si ninguna de estas distancias excede la tolerancia, la simplificación de vértices habrá terminado y quedarán en la simplificación los puntos inicial y final de la línea. • En el caso de que sí se exceda la tolerancia, el punto con mayor distancia se tomará como punto crítico que dividirá a la línea base en dos segmentos, y así, se seguirá repitiendo el proceso con cada segmento hasta que no sea necesario realizar más divisiones. A continuación, se muestra un diagrama de flujo del algoritmo de Ramer-Douglas-Peucker:
4.2 Problemas resueltos 63 Inicio: algoritmo de Ramer-Douglas-Peucker Polilínea a simplificar Tolerancia lineal Cálculo de puntos críticos a partir de la tolerancia lineal Cálculo de distancias perpendiculares a los puntos críticos ¿Excede alguna de las distancias perpendiculares la tolerancia lineal? División de la polilínea base en dos a partir del punto con mayor distancia a los puntos críticos Fin: algoritmo de Ramer-Douglas-Peucker Sí No Figura 4.3 Diagrama de flujo del algoritmo de Ramer-Douglas-Peucker. El algoritmo de Ramer-Douglas-Peucker se aplicará a la región tormentosa a través del comando reducem en MATLAB. Este, implementa el método recibiendo directamente los vértices en coordenadas sobre la Tierra (φ,λ) . Las tolerancias se escogerán ajustándolas al número de puntos que se pretenda conseguir. Del mismo modo, tras la aplicación del comando, es necesario comprobar que los polígonos resultantes no se solapen. En el caso de la región en consideración y para las tolerancias escogidas, esto no sucede. En conclusión, y pese a que el autor ha resuelto muchos más, se mostrarán seis problemas, cuya elección de puntos de origen y destino es variada, cada uno de ellos resuelto para cinco regiones con distinto número de puntos. Se mostrará el grafo de visibilidad únicamente para el primer problema resuelto debido a su similitud con el resto, así como la trayectoria más corta para todos los problemas en el caso del número de puntos original de la región. Tras lo anterior, se mostrarán los tiempos medios de ejecución resultantes de los seis problemas (de nuevo, no se muestran los tiempos para cada problema por separado debido a su similitud) y se realizará un análisis de los mismos. 4.2 Problemas resueltos Se adjuntan a continuación seis problemas resueltos para la puesta en marcha de los algoritmos y la comprobación de su correcto funcionamiento.
64 Capítulo 4. Resultados 4.2.1 Problema 1. Vuelo Bruselas-Sarajevo En este primer problema, el transcurso de la aeronave a lo largo de la región tormentosa se realizará de forma transversal. El grafo de visibilidad resultante se muestra a modo de ejemplo en la Figura 4.4: Figura 4.4 Grafo de visibilidad del vuelo Bruselas-Sarajevo. Asimismo, la trayectoria más corta obtenida se muestra en la Figura 4.5. Se tiene un camino desarrollado a lo largo de ocho nodos incluyendo origen y destino. Figura 4.5 Trayectoria más corta del vuelo Bruselas-Sarajevo.
4.2 Problemas resueltos 65 4.2.2 Problema 2. Vuelo Madrid-Kiev Se propone en este caso un problema en el que el transcurso de la aeronave a lo largo de la región tormentosa se realiza de forma longitudinal. La trayectoria más corta obtenida se muestra en la Figura 4.6. Se tiene un camino desarrollado a lo largo de diez nodos incluyendo origen y destino. Figura 4.6 Trayectoria más corta del vuelo Madrid-Kiev. 4.2.3 Problema 3. Vuelo Turín-Varsovia En esta situación, se tiene un problema cuyo punto de origen es una ciudad interior a la región tormentosa y cuyo punto de destino es una ciudad exterior. La trayectoria más corta obtenida se muestra en la Figura 4.7. Se tiene un camino desarrollado a lo largo de 14 nodos incluyendo origen y destino. Figura 4.7 Trayectoria más corta del vuelo Turín-Varsovia.
72 Capítulo 5. Conclusiones y líneas futuras podría contribuir a mejorar aún más la solución obtenida a través del método de Lee. • Implementación de algoritmos con un enfoque distinto al de la generación del grafo de visibilidad empleando todos los nodos de los polígonos obstáculo Tras el empleo de los algoritmos desarrollados en este trabajo, se obtiene un conjunto de conexiones entre nodos muy extenso. No obstante, un gran número de las trayectorias derivadas de estas conexiones podrían despreciarse debido a que, en situaciones reales, una aeronave nunca las seguiría. En consecuencia, sería interesante la implementación de algoritmos que generasen únicamente conexiones relevantes para la ejecución de los planes de vuelo de un avión y eliminen el costo de computación asociado a las trayectorias irrelevantes. El planteamiento de estos algoritmos se podría basar en el cálculo de las envolventes convexas de los polígonos que conforman la región tormentosa para, tras ello, calcular 4 tangentes para cada posible emparejamiento de polígonos, y, posteriormente, pasar a comprobar si son visibles o no (lo serán si no cortan un tercer polígono). Por otra parte, también se trazarían las tangentes desde el origen a todos los polígonos y del destino a todos los polígonos, procediéndose a eliminar aquellas que no sean visibles. Este planteamiento podría contribuir a disminuir significativamente el costo computacional del problema, eso sí, generando un número de conexiones entre nodos considerablemente menor.
Apéndice A Cálculo de distancias loxodrómicas En este apéndice, se muestra el procedimiento empleado en este trabajo para calcular las distancias loxodrómicas dlox entre los vértices mutuamente visibles del grafo. Con el fin de evitar el uso de comandos de MATLAB que conllevan un alto costo computacional, se va a recurrir a un planteamiento basado en una serie de fórmulas analíticas. Sean (φi,λi) y (φj,λj) las coordenadas de dos puntos i y j entre los que se quiere calcular la distancia loxodrómica, y sean RT el radio de la Tierra y h la altura de vuelo. El proceso para computar la distancia loxodrómica dlox es el siguiente ([5]): 1. Si φj=φi a) Cálculo del rumbo loxodrómico χlox: χlox =arctan λj−λi lntan(π/4−φi/2) tan(π/4−φj/2) Si λj−λi≤π(A.1) χlox =arctan λj+2π−λi lntan(π/4−φi/2) tan(π/4−φj/2) Si λj−λi>πyλj<0(A.2) χlox =arctan λj−2π−λi lntan(π/4−φi/2) tan(π/4−φj/2) Si λj−λi>πyλj>0(A.3) •Corrección del rumbo loxodrómico: χlox =π+χlox Si φj<φi(A.4) b) Cálculo de la distancia angular loxodrómica αlox: αlox =φj−φi cos χlox (A.5) 73
74 Apéndice A. Cálculo de distancias loxodrómicas 2. Si φi=φj=φ αlox =λj−λicosφ(A.6) 3. Cálculo de la distancia loxodrómica: dlox =αlox(RT+h)(A.7) Nótese que, en este trabajo, puesto que se considera una altitud de vuelo constante para la generación de la trayectoria más corta, se ha tomado h=0, pues no es influyente en el resultado. A continuación se muestra la función implementada en MATLAB para el cálculo de las distancias loxodrómicas. 1function d_lox = LoxodromicDistance(phi_i, lambda_i, phi_j, lambda_j, h) 2 3R_T = 6378.14; % Radio de la Tierra en kilómetros 4 5if phi_j == phi_i 6alpha_lox = abs(lambda_j - lambda_i)*cos(phi_i); 7 8else 9if abs(lambda_j - lambda_i) > pi 10 if lambda_j < 0 11 ji_lox = atan((lambda_j + 2*pi - lambda_i)/log((tan(pi/4 - phi_i/2)/... 12 tan(pi/4 - phi_j/2)))); 13 else 14 ji_lox = atan((lambda_j - 2*pi - lambda_i)/log((tan(pi/4 - phi_i/2)/... 15 tan(pi/4 - phi_j/2)))); 16 end 17 else 18 ji_lox = atan((lambda_j - lambda_i)/log((tan(pi/4 - phi_i/2)/tan(pi/4 - ... 19 phi_j/2)))); 20 end 21 22 if phi_j < phi_i 23 ji_lox = ji_lox + pi; 24 end 25 26 alpha_lox = (phi_j - phi_i)/cos(ji_lox); 27 28 end 29 30 d_lox = (R_T + h)*alpha_lox; 31 32 end Código A.1 Función para el cálculo de distancias loxodrómicas.
Apéndice B Algoritmos de ordenación En este apéndice, se muestran los algoritmos de ordenación empleados cuando ha sido requerido en el desarrollo del método de Lee. Estos, se basan en el método HeapSort que garantiza que, para realizar una ordenación de n elementos, el costo computacional de la misma en el peor de los casos es de O(nlogn). Puesto que el objetivo de este trabajo no es la codificación del método, HeapSort ([ 6 ]), no se va a explicar su funcionamiento. No obstante, sí se va a señalar que los algoritmos HeapSort usualmente programados ordenan un vector según el valor de sus elementos, mientras que los métodos HeapSort empleados en este trabajo han empleado otras condiciones para la ordenación de sus elementos. Estas se detallan en las siguientes secciones. B.1 Primera versión del método HeapSort La primera versión del algoritmo de ordenación HeapSort se ha empleado para realizar la ordenación de vértices en sentido antihorario a partir del sentido positivo del eje x alrededor de un nodo i . Para ello, se ha hecho uso de las coordenadas de los vértices alrededor del vértice i(citando [6]). Considérese dos nodos genéricos w1 y w2 tales que w1 se encuentra en una posición inferior que w2 en la lista donde se han de ordenar los vértices de los distintos nodos alrededor del nodo i . El procedimiento seguido para decidir si se han de permutar las posiciones de w1 y w2 se muestra en el Pseudocódigo B.1. En primer lugar, se hacen comprobaciones basadas en el incremento de la coordenada y de los puntos w1 y w2 con respecto a la coordenada y del nodo i . Por otra parte, se comparan los nodos w1 y w2 cuando tienen la misma coordenada y que el vértice i . Si lo anterior no permite la ordenación de los nodos, se comprueban las desigualdades relativas al producto vectorial de la componente no nula de lo que serían los vectores iw1 y iw2 . En caso de que estos vectores tengan la misma dirección, se comprueba la distancia de w1yw2al nodo i. 75
76 Capítulo B. Algoritmos de ordenación Pseudocódigo B.1 Less(w1,w2,i) Entrada: Las coordenadas de los nodos entre los que se plantea la permutación, (w1,x,w1,y) y (w2,x,w2,y) , y las coordenadas del nodo alrededor del que se realiza la ordenación de vértices, (ix,iy). Salida: Verdadero oFalso. 1: si w2,y−iy≥0yw1,y−iy<0entonces 2: devolver verdadero 3: si no si w2,y−iy<0yw1,y−iy≥0entonces 4: devolver falso 5: si no si w2,y−iy=0yw1,y−iy=0entonces 6: si w2,x−ix>0entonces 7: si w1,x−ix>0entonces 8: si w2,x<w1,xentonces 9: devolver verdadero 10: si no 11: devolver falso 12: si no 13: devolver verdadero 14: si no 15: si w1,x−ix<0entonces 16: si w2,x>w1,xentonces 17: devolver verdadero 18: si no 19: devolver falso 20: si no 21: devolver falso 22: si no 23: det = (w2,x−ix)(w1,y−iy)−(w2,y−iy)(w1,x−ix) 24: si det >0entonces 25: devolver verdadero 26: si no si det <0entonces 27: devolver falso 28: si no 29: si (w2,x−ix)2+(w2,y−iy)2<(w1,x−ix)2+(w1,y−iy)2entonces 30: devolver verdadero 31: si no 32: devolver falso La primera versión de la función HeapSort recibirá como argumentos de entrada una matriz A (nótese que en este caso, A no es uno de los nodos extremos de e1 como se mostró en la sección 3.2.3) que se corresponderá con la matriz mostrada en 3.7 y, además, las coordenadas del nodo i alrededor del que se realiza la ordenación, en este caso denotado como c. A continuación, se muestra el código de las funciones desarrolladas para la implementación de la primera versión del método HeapSort en este trabajo. 1function A = HeapSort_v1(A, c) 2Dim = size(A); N = Dim(1);
B.1 Primera versión del método HeapSort 77 3LastTopPos = floor(N/2); 4 5for i = LastTopPos:-1:1 6A = Heapify_v1(A, N, i, c); 7end 8 9for i = N-1:-1:1 10 aux = A(i+1,:); 11 A(i+1,:) = A(1,:); 12 A(1,:) = aux; 13 A = Heapify_v1(A, i, 1, c); 14 end 15 end Código B.1 Primera versión del algoritmo HeapSort. 1function A = Heapify_v1(A, N, i, c) 2BiggerPos = i; 3LeftPos = 2*i; 4RightPos = 2*i + 1; 5 6if LeftPos <= N 7if A(BiggerPos,2) >= c(2) && A(LeftPos,2) < c(2) 8flag = 1; 9 10 elseif A(BiggerPos,2) < c(2) && A(LeftPos,2) >= c(2) 11 flag = 0; 12 13 elseif A(BiggerPos,2) == c(2) && A(LeftPos,2) == c(2) 14 if A(BiggerPos,1) > c(1) 15 if A(LeftPos,1) > c(1) 16 if A(BiggerPos,1) < A(LeftPos,1) 17 flag = 1; 18 else 19 flag = 0; 20 end 21 else 22 flag = 1; 23 end 24 25 else 26 if A(LeftPos,1) < c(1) 27 if A(BiggerPos,1) > A(LeftPos,1) 28 flag = 1; 29 else 30 flag = 0; 31 end 32 else 33 flag = 0; 34 end 35 end 36 37 else 38 det = (A(BiggerPos,1) - c(1))*(A(LeftPos,2) - c(2)) - (A(BiggerPos,2) - ... 39 c(2))*(A(LeftPos,1) - c(1)); 40 41 if det > 0 42 flag = 1; 43 elseif det < 0
78 Capítulo B. Algoritmos de ordenación 44 flag = 0; 45 else 46 if (A(BiggerPos,1) - c(1))^2 + (A(BiggerPos,2) - c(2))^2 ... 47 < (A(LeftPos,1) - c(1))^2 + (A(LeftPos,2) - c(2))^2 48 flag = 1; 49 else 50 flag = 0; 51 end 52 end 53 end 54 if flag == 1 55 BiggerPos = LeftPos; 56 end 57 end 58 59 if RightPos <= N 60 if A(BiggerPos,2) >= c(2) && A(RightPos,2) < c(2) 61 flag = 1; 62 63 elseif A(BiggerPos,2) < c(2) && A(RightPos,2) >= c(2) 64 flag = 0; 65 66 elseif A(BiggerPos,2) == c(2) && A(RightPos,2) == c(2) 67 if A(BiggerPos,1) > c(1) 68 if A(RightPos,1) > c(1) 69 if A(BiggerPos,1) < A(RightPos,1) 70 flag = 1; 71 else 72 flag = 0; 73 end 74 else 75 flag = 1; 76 end 77 78 else 79 if A(RightPos,1) < c(1) 80 if A(BiggerPos,1) > A(RightPos,1) 81 flag = 1; 82 else 83 flag = 0; 84 end 85 else 86 flag = 0; 87 end 88 end 89 90 else 91 det = (A(BiggerPos,1) - c(1))*(A(RightPos,2) - c(2)) - (A(BiggerPos,2) - ... 92 c(2))*(A(RightPos,1) - c(1)); 93 94 if det > 0 95 flag = 1; 96 elseif det < 0 97 flag = 0; 98 else 99 if (A(BiggerPos,1) - c(1))^2 + (A(BiggerPos,2) - c(2))^2 < ... 100 (A(RightPos,1) - c(1))^2 + (A(RightPos,2) - c(2))^2 101 flag = 1; 102 else 103 flag = 0; 104 end 105 end
B.2 Segunda versión del método HeapSort 79 106 end 107 if flag == 1 108 BiggerPos = RightPos; 109 end 110 end 111 112 if BiggerPos ~= i 113 aux = A(i,:); 114 A(i,:) = A(BiggerPos,:); 115 A(BiggerPos,:) = aux; 116 117 A = Heapify_v1(A, N, BiggerPos, c); 118 end 119 end Código B.2 Primera versión de la función auxiliar Heapify del algoritmo HeapSort. B.2 Segunda versión del método HeapSort La segunda versión del algoritmo de ordenación HeapSort se ha empleado para realizar la ordenación de las aristas que se encuentran almacenadas en el conjunto ordenado T . Para ello, se ha hecho uso tanto de las distancias de las aristas al nodo i almacenadas en ρ como del conjunto xcomp que almacena la componente en el eje x del vector característico normalizado de las aristas para desempatar en cuanto a la ordenación en el caso de que dos aristas presentasen la misma distancia al nodo i . En este caso, las condiciones para la ordenación son sencillas: 1. Las aristas con menor distancia al nodo iirán primero. 2. En caso de empates, la arista con menor componente en el eje x de su vector característico normalizado irá primero. La segunda versión de la función HeapSort recibirá como argumento de entrada una matriz A que, en este caso, se corresponde con la mostrada en 3.9. A continuación, se muestra el código de las funciones desarrolladas para la implementación de la segunda versión del método HeapSort en este trabajo. 1function A = HeapSort_v2(A) 2Dim = size(A); N = Dim(1); 3LastTopPos = floor(N/2); 4 5for i = LastTopPos:-1:1 6A = Heapify_v2(A, N, i); 7end 8 9for i = N-1:-1:1 10 aux = A(i+1,:); 11 A(i+1,:) = A(1,:); 12 A(1,:) = aux; 13 A = Heapify_v2(A, i, 1); 14 end
80 Capítulo B. Algoritmos de ordenación 15 end Código B.3 Segunda versión del algoritmo HeapSort. 1function A = Heapify_v2(A, N, i) 2rho = A(:, length(A(1,:))); 3x_comp = A(:, length(A(1,:)) - 1); 4BiggerPos = i; 5LeftPos = 2*i; 6RightPos = 2*i + 1; 7 8if [LeftPos <= N && rho(BiggerPos) < rho(LeftPos)] || [LeftPos <= N && ... 9rho(BiggerPos) == rho(LeftPos) && x_comp(BiggerPos) < x_comp(LeftPos)] 10 BiggerPos = LeftPos; 11 end 12 13 if [RightPos <= N && rho(BiggerPos) < rho(RightPos)] || [RightPos <=N && ... 14 rho(BiggerPos) == rho(RightPos) && x_comp(BiggerPos) < x_comp(RightPos)] 15 BiggerPos = RightPos; 16 end 17 18 if BiggerPos ~= i 19 aux = A(i,:); 20 A(i,:) = A(BiggerPos,:); 21 A(BiggerPos,:) = aux; 22 23 A = Heapify_v2(A, N, BiggerPos); 24 end 25 end Código B.4 Segunda versión de la de la función auxiliar Heapify del algoritmo HeapSort.
Apéndice C Algoritmo de Búsqueda Binaria e Inserción En este apéndice, se muestra el algoritmo responsable de realizar la correcta inserción del conjunto de aristas obstáculo dadas en el conjunto newedges en T cuando ha sido requerido en el desarrollo del método de Lee con un costo computacional, en el peor de los casos, de O(log2nT) . Puesto que el objetivo de este trabajo no es la codificación del algoritmo de Búsqueda Binaria, no se va a explicar su implementación. No obstante, sí se va a señalar la diferencia del algoritmo empleado con respecto a los métodos de Búsqueda Binaria convencionales ([7]). Un algoritmo de Búsqueda Binaria usual encuentra la posición de un determinado valor en un vector ordenado. La necesidad de este trabajo es encontrar en un determinado conjunto ( ρ ) la posición a insertar newedges en T, así como su distancia al nodo ien ρ. Nótese que a lo largo de los algoritmos VisibleVerticesVA (i,V) yVisibleVerticesOD (i,V) nunca se actualiza la distancia de las aristas obstáculo almacenadas en T al nodo i desde su inserción inicial. Debido a ello, en este algoritmo será necesario recalcular la distancia de las aristas que se encuentren introducidas en T al nodo i , pero únicamente de aquellas con las que el algoritmo vaya realizando comprobaciones. Para actualizar esta distancia, se volverá a aplicar el procedimiento descrito en la sección 3.2.3 basado en la Figura 3.15, donde, en este caso, e1 sería la arista de T con la que se está realizando la comparación, y la búsqueda de las coordenadas de los nodos A y B se realizará directamente en la matriz VMercator. A continuación, se muestra el código de la función desarrollada para la implementación del algoritmo de Búsqueda Binaria e Inserción en este trabajo. Nótese que en este algoritmo se emplean índices i y j pero su interpretación no tiene nada que ver con aquella asociada a los índices de los nodos que se ha estado empleando a lo largo del trabajo. 1function [Pos, rho, T] = BinarySearchAndInsertion(rho, rho_k0, T, newedges, ... 2V_Mercator, origin, destiny) 3n_T = length(rho); 4i = 1; 5j = n_T; 6while i <= j 7medium = fix((i+j)/2); 81
88 Capítulo D. Códigos de MATLAB 17 % del número de polígonos obstáculo del grafo 18 W=find(isnan(V(:,1))); 19 npol = length(W) - 1; 20 21 % Matriz de WPs, orden del grafo y número de aristas obstáculo del grafo 22 V_0 = V; V_0(W,:) = []; % Primero se suprimen los NaN 23 WPs = unique(V_0, ‘rows’, ‘stable’); % Seguidamente se suprimen los vértices repetidos 24 N=length(WPs); n = N - 2; 25 26 % Implementación de la transformación de Mercator 27 Lon_0 = mean(WPs(:,2)); 28 x_V = (V(:,2) - Lon_0); y_V = log(tan(pi/4 + V(:,1)/2)); V_Mercator = [x_V, y_V]; 29 x_WPs = (WPs(:,2) - Lon_0); y_WPs = log(tan(pi/4 + WPs(:,1)/2)); WPs_Mercator = ... 30 [x_WPs, y_WPs]; 31 32 % Vector de índices de los puntos del grafo 33 index = 1:N; 34 35 % Cálculo de las matrices E y Einv 36 E=NaN(length(V),2); Einv = zeros(N,2); 37 for p = 1:npol 38 if p == 1 39 E(1: W(p) - 2,:) = [1 : W(p) - 2; 2 : W(p) - 2, 1]’; 40 Einv(1: W(p) - 2,:) = [1 : W(p) - 2; W(p) - 2, 1 : W(p) - 3]’; 41 else 42 E(W(p-1) + 1: W(p) - 2,:) = [W(p-1) + 1 - 2*(p-1) : W(p) - 2*p; ... 43 W(p-1) + 2 - 2*(p-1) : W(p) - 2*p, W(p-1) + 1 - 2*(p-1)]’; 44 Einv(W(p-1) - 2*(p-1) + 1 : W(p) - 2*p,:) = [W(p-1) + 1 : W(p) - 2; ... 45 W(p) - 2, W(p-1) + 1 : W(p) - 3]’; 46 end 47 end 48 49 % Inicialización de la matriz de adyacencia lógica 50 Conex = false(N,N); 51 52 %% CÓMPUTO DEL GRAFO DE VISIBILIDAD 53 for i = 1:N 54 if i <= n 55 Conex = VisibleVerticesVA(i, V_Mercator(1:N_in - 4,:), WPs_Mercator(1:n,:), ... 56 index(1:n), E, Einv, Conex, N_in); 57 else 58 Conex = VisibleVerticesOD(i, V_Mercator, WPs_Mercator, N, index, E, Einv, ... 59 Conex, N_in); 60 end 61 end 62 63 %% OBTENCIÓN DEL CAMINO MÁS CORTO ENTRE O Y D 64 65 % Creación de la matriz de adyacencia con pesos 66 [i_vec, j_vec] = find(Conex); 67 for n_c = 1:length(i_vec) 68 d_lox(n_c) = LoxodromicDistance(WPs(i_vec(n_c),1), WPs(i_vec(n_c),2), ... 69 WPs(j_vec(n_c),1), WPs(j_vec(n_c),2), 0); 70 end 71 72 significant_numbers = 9; % Número de cifras significativas que se considera 73 factor = 10.^(significant_numbers - ceil(log10(abs(d_lox)))); 74 d_lox = round(d_lox .* factor) ./ factor; 75 76 Adj = sparse(i_vec,j_vec,d_lox); 77 78 % Cálculo del camino más corto
D.3 Algoritmo de Lee 89 79 G = graph(Adj); 80 [path, distance] = shortestpath(G, N - 1, N); 81 82 %% OPERACIONES AUXILIARES AL ALGORITMO 83 84 %%% Para trazar el grafo de visibilidad, descomentar las siguientes líneas de código 85 % figure(1) 86 %fori=1:N 87 % for j = i : N 88 % if Conex(i,j) == 1 89 % geoplot([WPs(i,1), WPs(j,1)]*180/pi, [WPs(i,2), WPs(j,2)]*180/pi); 90 % hold on 91 % end 92 % end 93 % end 94 % 95 %%% Para trazar el camino más corto, descomentar las siguientes líneas de código 96 % figure(2) 97 % geoplot(V(:,1)*180/pi, V(:,2)*180/pi,’b’); 98 % hold on 99 % geoplot(WPs(path,1)*180/pi, WPs(path,2)*180/pi,’r’); 100 % hold on 101 % for i=1:length(path) 102 % geoplot(WPs(path(i),1)*180/pi,WPs(path(i),2)*180/pi,’rd’) 103 % end Código D.3 Programa principal del algoritmo de Lee. D.3.2 Función VISIBLEVERTICESVA(i,V) 1function Conex = VisibleVerticesVA(i, V_Mercator, WPs_Mercator, N, index, E, Einv, ... 2Conex, N_in) 3 4%% ORDENACIÓN DE VÉRTICES EN SENTIDO ANTIHORARIO ALREDEDOR DEL NODO i 5 6VerticesData = [WPs_Mercator, index’]; 7VerticesData(i,:) = []; 8VerticesData((VerticesData(:,2) - WPs_Mercator(i,2)) < 0,:) = []; 9% 10 % Caso en el que solo queda un vértice restante, para el que se comprobará directamente 11 % la visibilidad 12 if length(VerticesData(:,3)) == 1 13 14 j = VerticesData(3); 15 16 % Cálculo del vector v 17 v = (WPs_Mercator(j, :) - WPs_Mercator(i, :))’; 18 19 % Formación de la base no ortogonal de vectores 20 r = (WPs_Mercator(E(Einv(i,1), 2),:) - WPs_Mercator(i,:))’; 21 s = (WPs_Mercator(E(Einv(i,2), 1),:) - WPs_Mercator(i,:))’; 22 23 % Producto vectorial de r y s 24 q_z = r(1)*s(2) - r(2)*s(1); 25 26 % Resolución del sistema de ecuaciones 27 a = (v(1)*s(2) - v(2)*s(1))/q_z; 28 b = (r(1)*v(2) - r(2)*v(1))/q_z;
90 Capítulo D. Códigos de MATLAB 29 30 % Aplicación del criterio con una cierta tolerancia 31 if (sign(q_z) == -1 && ~(a>1e-12 && b>1e-12)] || (sign(q_z) == 1 && ... 32 (a>=-1e-12 && b>=-1e-12)) 33 Conex(i,j) = 1; Conex(j,i) = 1; 34 end 35 36 % Caso en que queda más de un vértice, para el que se recorrerá el algoritmo completo 37 elseif length(VerticesData(:,3)) > 1 38 39 %% ORDENACIÓN DE VÉRTICES EN SENTIDO ANTIHORARIO 40 VerticesData= HeapSort_v1(VerticesData, WPs_Mercator(i,:)); 41 j = VerticesData(:,3); 42 43 %% INICIALIZACIÓN DEL CONJUNTO ORDENADO T 44 45 % Intersección del segmento sigmagorro con todos los polígonos 46 [xk, yk, K] = polyxpoly([WPs_Mercator(i,1), WPs_Mercator(i,1) + 2*pi], ... 47 [WPs_Mercator(i,2), WPs_Mercator(i,2)], V_Mercator(:,1), V_Mercator(:,2)); 48 49 % Inicialización de estructuras para introducir el conjunto ordenado T, distancia 50 % de las aristas en T al punto i en rho, posición de las aristas de T en Tinv y 51 % la componente x del vector característico de las aristas en x_comp 52 T=zeros(0,1); rho = zeros(0,1); Tinv = NaN(N_in,1); x_comp = zeros(0,1); 53 54 % Bucle para tratar las aristas que se deban insertar en T 55 for k0=1:length(xk) 56 57 % Primera comprobación: que el punto de intersección no sea el punto i del que 58 % parte el segmento sigmagorro 59 if xk(k0) ~= WPs_Mercator(i,1) 60 61 % Arista de intersección 62 k = K(k0,2); 63 64 % Segunda comprobación: si la intersección es justamente un vértice o si es 65 % un vértice repetido 66 isvertice = 0; % Variable que se igualará a 1 si se detecta el caso de 67 % vértice 68 isverticerep = 0; % Variable que se igualará a 1 si se detecta el caso de 69 % vértice repetido 70 71 if V_Mercator(k,1) == xk(k0) && V_Mercator(k,2) == yk(k0) 72 isvertice = 1; 73 elseif V_Mercator(k+1,1) == xk(k0) && V_Mercator(k+1,2) == yk(k0) 74 isverticerep = 1; 75 end 76 77 if isvertice == 1 78 % Caso en que la intersección es un vértice 79 kminus = Einv(E(k,1),2); 80 81 % Cálculo del vector t 82 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 83 84 % Para incluir la arista o no en T, se observa la componente en el eje 85 % y del vector t 86 if t_0(2) > 0 87 88 % Introducción de la arista correspondiente en T 89 T = [T; kminus]; 90
D.3 Algoritmo de Lee 91 91 % Cálculo de la distancia de la arista al punto i 92 rho_k0 = xk(k0) - WPs_Mercator(i,1); 93 rho = [rho; rho_k0]; 94 95 % Almacenamiento de la componente x del vector t normalizado por si 96 % se da el caso de que haya que desempatar en cuanto a la 97 % ordenación en T 98 x_comp = [x_comp; t_0(1)]; 99 end 100 101 % Las operaciones con el vector u son análogas a las realizadas con el 102 % vector t 103 u = V_Mercator(k + 1,:) - V_Mercator(k,:); u_0 = u/norm(u); 104 if u_0(2) > 0 105 T = [T; k]; 106 rho_k0 = xk(k0) - WPs_Mercator(i,1); 107 rho = [rho; rho_k0]; 108 x_comp = [x_comp; u_0(1)]; 109 end 110 111 elseif isverticerep ~= 1 112 % La intersección ya no puede ser un vértice 113 T = [T; k]; 114 rho_k0 = xk(k0) - WPs_Mercator(i,1); 115 rho = [rho; rho_k0]; 116 x_comp = [x_comp; NaN]; % Puesto que en este caso no será necesario 117 % desempatar en cuanto a la ordenación en T, se introduce un NaN 118 end 119 end 120 end 121 122 TreeData = [T, x_comp, rho]; 123 124 % Ordenación a través del método HeapSort 125 TreeData = HeapSort_v2(TreeData); 126 T = TreeData(:,1); rho = TreeData(:,3); 127 128 % Relleno del conjunto Tinv 129 for y = 1:length(T) 130 Tinv(T(y)) = y; 131 end 132 133 %% COMPROBACIÓN DE LA VISIBILIDAD DE LOS VÉRTICES 134 135 for j0 = 1:length(j) 136 if j0 == 1 || (j0 > 1 && (WPs_Mercator(j(j0-1),1) - WPs_Mercator(i,1))* ... 137 (WPs_Mercator(j(j0),2) - WPs_Mercator(i,2)) ~= ... 138 (WPs_Mercator(j(j0),1) - WPs_Mercator(i,1))* ... 139 (WPs_Mercator(j(j0-1),2) - WPs_Mercator(i,2))) % No se considerarán 140 % casos redundantes 141 142 if isempty(T) 143 vis = 1; 144 else 145 vis = VisibleVA(j(j0), i, WPs_Mercator, T(1), E, Einv); 146 end 147 if vis == 1 148 Conex(i,j(j0)) = 1; Conex(j(j0),i) = 1; 149 end 150 end 151 152 if ~(WPs_Mercator(j(j0), 2) == WPs_Mercator(i,2) && WPs_Mercator(j(j0),1) ...
92 Capítulo D. Códigos de MATLAB 153 > WPs_Mercator(i,1)) && j0 < length(j) % Porque en el primer caso ya se 154 % ha realizado la inicialización de T y para el último vértice no es 155 % necesario actualizar % T pues no se comprobará la visibilidad 156 % de ningún otro nodo 157 158 %% ACTUALIZACIÓN DE T 159 160 % Búsqueda de la existencia en T tanto de la arista k como de la arista 161 % k^{-} 162 k = Einv(j(j0),1); kminus = Einv(j(j0),2); 163 pos_k = Tinv(k); pos_kminus = Tinv(kminus); 164 165 % Caso en que sea necesario eliminar ambas aristas en T 166 if ~isnan(pos_k) && ~isnan(pos_kminus) 167 168 Tinv([k; kminus]) = NaN; 169 if pos_kminus > pos_k 170 Tinv(T(pos_kminus + 1:length(T))) = Tinv(T(pos_kminus + ... 171 1:length(T))) - 2; 172 else 173 Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1:length(T))) - 2; 174 end 175 T([pos_k; pos_kminus]) = []; 176 rho([pos_k; pos_kminus]) = []; 177 178 % Casos en que sea necesario eliminar una arista e introducir otra (si esta 179 % no va en la dirección del segmento ij_{j0}) 180 elseif ~isnan(pos_k) && isnan(pos_kminus) 181 182 % Vector director del segmento ij_{j0} 183 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 184 185 % Vector t 186 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 187 188 if (abs(t_0(1) - v_0(1)) < 1e-12 && abs(t_0(2) - v_0(2)) < 1e-12) ... 189 || (abs(t_0(1) + v_0(1)) < 1e-12 && abs(t_0(2) + v_0(2)) ... 190 < 1e-12) % Se ha aplicado una cierta tolerancia 191 % Caso en que solo se suprime una arista 192 Tinv(k) = NaN; Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1:... 193 length(T))) - 1; 194 T(pos_k) = []; 195 rho(pos_k) = []; 196 else 197 % Caso en que se sustituye una arista por la otra 198 Tinv(k) = NaN; Tinv(kminus) = pos_k; 199 T(pos_k) = kminus; 200 end 201 202 % Caso análogo al anterior 203 elseif isnan(pos_k) && ~isnan(pos_kminus) 204 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 205 u = V_Mercator(k+1,:) - V_Mercator(k,:); u_0 = u/norm(u); 206 207 if (abs(u_0(1) - v_0(1)) < 1e-12 && abs(u_0(2) - v_0(2)) < 1e-12) ... 208 || (abs(u_0(1) + v_0(1)) < 1e-12 && abs(u_0(2) + v_0(2)) ... 209 < 1e-12) 210 Tinv(kminus) = NaN; Tinv(T(pos_kminus + 1:length(T))) = ... 211 Tinv(T(pos_kminus + 1:length(T))) - 1; 212 T(pos_kminus) = []; 213 rho(pos_kminus) = []; 214 else
D.3 Algoritmo de Lee 93 215 Tinv(kminus) = NaN; Tinv(k) = pos_kminus; 216 T(pos_kminus) = k; 217 end 218 219 % Caso en que hay que añadir dos aristas (si la dirección de estas no va 220 % en la del segmento ij_{j0}) 221 else 222 223 % Vector director del segmento ij_{j0} 224 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 225 226 % Vectores u y t 227 u = V_Mercator(k + 1,:) - V_Mercator(k,:); u_0 = u/norm(u); 228 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 229 230 % Distancia al nodo i 231 rho_k0 = norm(v); 232 233 % Definición de las aristas a introducir en T 234 if (abs(u_0(1) - v_0(1)) < 1e-12 && abs(u_0(2) - v_0(2)) < 1e-12) ... 235 || (abs(u_0(1) + v_0(1)) < 1e-12 && abs(u_0(2) + v_0(2)) ... 236 < 1e-12) 237 newedges = kminus; % Solo se añadirá una arista 238 239 elseif (abs(t_0(1) - v_0(1)) < 1e-12 && abs(t_0(2) - v_0(2)) ... 240 < 1e-14) || (abs(t_0(1) + v_0(1)) < 1e-14 && abs(t_0(2) + ... 241 v_0(2)) < 1e-12) 242 newedges = k; % Solo se añadirá una arista 243 244 else 245 246 %Se añadirán dos aristas 247 rho_k0 = [rho_k0; rho_k0]; 248 249 % Las dos aristas se ordenan a través del producto escalar 250 udotv = dot(u_0,v); 251 tdotv = dot(t_0,v); 252 253 if udotv > tdotv 254 newedges = [kminus; k]; 255 else 256 newedges = [k; kminus]; 257 end 258 259 end 260 261 % Introducción en T y actualización de rho y Tinv 262 if isempty(T) 263 T = newedges; 264 rho = rho_k0; 265 Tinv(T(1:length(newedges))) = 1:length(newedges); 266 else 267 [pos_newedges, rho, T] = BinarySearchAndInsertion(rho, rho_k0, ... 268 T, newedges, V_Mercator, WPs_Mercator(i,:), ... 269 WPs_Mercator(j(j0),:)); 270 Tinv(T(pos_newedges)) = pos_newedges; 271 Tinv(T(pos_newedges(length(pos_newedges)) + 1:length(T))) = ... 272 Tinv(T(pos_newedges(length(pos_newedges)) + 1:length(T))) ... 273 +length(pos_newedges); 274 end 275 end 276 end
94 Capítulo D. Códigos de MATLAB 277 end 278 end 279 end Código D.4 Función VisibleVerticesVA(i,V). D.3.3 Función VISIBLEVERTICESOD(i,V) 1function Conex = VisibleVerticesOD(i, V_Mercator, WPs_Mercator, N, index, E, Einv, ... 2Conex, N_in) 3 4%% ORDENACIÓN DE VÉRTICES EN SENTIDO ANTIHORARIO 5 6VerticesData = [WPs_Mercator, index’]; 7VerticesData(i,:) = []; 8 9% En esta función, siempre será necesario recorrer el algoritmo completo 10 VerticesData= HeapSort_v1(VerticesData, WPs_Mercator(i,:)); 11 j = VerticesData(:,3); 12 13 %% INICIALIZACIÓN DEL CONJUNTO ORDENADO T 14 15 % Intersección del segmento sigmagorro con todos los polígonos 16 [xk, yk, K] = polyxpoly([WPs_Mercator(i,1), WPs_Mercator(i,1) + 2*pi], ... 17 [WPs_Mercator(i,2), WPs_Mercator(i,2)], V_Mercator(:,1), V_Mercator(:,2)); 18 19 % Inicialización de estructuras para introducir el conjunto ordenado T, distancia de 20 % las aristas en T al punto i en rho, posición de las aristas de T en Tinv y la 21 % componente x del vector característico de las aristas en x_comp 22 T=zeros(0,1); rho = zeros(0,1); Tinv = NaN(N_in,1); x_comp = zeros(0,1); 23 24 % Bucle para tratar las aristas que se deban insertar en T 25 for k0=1:length(xk) 26 27 N=length(WPs_Mercator); 28 % En este caso la primera comprobación es que el punto de intersección no ni sea el 29 % nodo de origen O ni el de destino D 30 if xk(k0) ~= WPs_Mercator(N, 1) && xk(k0) ~= WPs_Mercator(N-1, 1) 31 32 % Arista de intersección 33 k = K(k0,2); 34 35 % Segunda comprobación: si la intersección es justamente un vértice o si es un 36 % vértice repetido 37 isvertice = 0; % Variable que se igualará a 1 si se detecta el caso de vértice 38 isverticerep = 0; % Variable que se igualará a 1 si se detecta el caso de 39 % vértice repetido 40 41 if V_Mercator(k,1) == xk(k0) && V_Mercator(k,2) == yk(k0) 42 isvertice = 1; 43 elseif V_Mercator(k + 1,1) == xk(k0) && V_Mercator(k + 1,2) == yk(k0) 44 isverticerep = 1; 45 end 46 47 if isvertice == 1 48 % Caso en que la intersección es un vértice 49 kminus = Einv(E(k,1),2); 50
D.3 Algoritmo de Lee 95 51 % Cálculo del vector t 52 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 53 54 % Para incluir la arista o no en T, se observa la componente en el eje y 55 % del vector t 56 if t_0(2) > 0 57 58 % Introducción de la arista correspondiente en T 59 T = [T; kminus]; 60 61 % Cálculo de la distancia de la arista al punto i 62 rho_k0 = xk(k0) - WPs_Mercator(i,1); 63 rho = [rho; rho_k0]; 64 65 % Almacenamiento de la componente x del vector t normalizado por si se 66 % da el caso de que haya que desempatar en cuanto a la ordenación en T 67 x_comp = [x_comp; t_0(1)]; 68 end 69 70 % Las operaciones con el vector u son análogas a las realizadas con el 71 % vector t 72 u = V_Mercator(k + 1,:) - V_Mercator(k,:); u_0 = u/norm(u); 73 if u_0(2) > 0 74 T = [T; k]; 75 rho_k0 = xk(k0) - WPs_Mercator(i,1); 76 rho = [rho; rho_k0]; 77 x_comp = [x_comp; u_0(1)]; 78 end 79 80 elseif isverticerep ~= 1 81 % La intersección ya no puede ser un vértice 82 T = [T; k]; 83 rho_k0 = xk(k0) - WPs_Mercator(i,1); 84 rho = [rho; rho_k0]; 85 x_comp = [x_comp; NaN]; % Puesto que en este caso no será necesario 86 % desempatar en cuanto a la ordenación en T, se introduce un NaN 87 end 88 end 89 end 90 91 TreeData = [T, x_comp, rho]; 92 93 % Ordenación a través del método HeapSort 94 TreeData = HeapSort_v2(TreeData); 95 T = TreeData(:,1); rho = TreeData(:,3); 96 97 % Se rellena el conjunto Tinv 98 for y = 1:length(T) 99 Tinv(T(y)) = y; 100 end 101 102 %% COMPROBACIÓN DE LA VISIBILIDAD DE LOS VÉRTICES 103 104 for j0 = 1:length(j) 105 if j0 == 1 || (j0 > 1 && (WPs_Mercator(j(j0-1),1) - WPs_Mercator(i,1))*... 106 (WPs_Mercator(j(j0),2) - WPs_Mercator(i,2)) ~= (WPs_Mercator(j(j0),1) - ... 107 WPs_Mercator(i,1))*(WPs_Mercator(j(j0-1),2) - WPs_Mercator(i,2))) % No se 108 % considerarán casos redundantes 109 110 if isempty(T) 111 vis = 1; 112 else
96 Capítulo D. Códigos de MATLAB 113 vis = VisibleOD(j(j0), i, WPs_Mercator, T(1), E); 114 end 115 end 116 if vis == 1 117 Conex(i,j(j0)) = 1; Conex(j(j0),i) = 1; 118 end 119 120 if (WPs_Mercator(j(j0), 2) ~= WPs_Mercator(i,2) || WPs_Mercator(j(j0),1) < ... 121 WPs_Mercator(i,1)) && j0 < length(j) && j(j0) ~= N && j(j0) ~= N-1 % Porque 122 % en el primer caso ya se ha realizado la inicialización de T y para el 123 % último vértice no es necesario actualizar T pues no se comprobará la 124 % visibilidad de ningún otro vértice 125 126 %% ACTUALIZACIÓN DE T 127 128 % Búsqueda de la existencia en T tanto de la arista k como de la arista k^{-} 129 k = Einv(j(j0),1); kminus = Einv(j(j0),2); 130 pos_k = Tinv(k); pos_kminus = Tinv(kminus); 131 132 % Caso en que sea necesario eliminar ambas aristas en T 133 if ~isnan(pos_k) && ~isnan(pos_kminus) 134 135 Tinv([k; kminus]) = NaN; 136 if pos_kminus > pos_k 137 Tinv(T(pos_kminus + 1:length(T))) = Tinv(T(pos_kminus + 1: ... 138 length(T))) - 2; 139 else 140 Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1:length(T))) - 2; 141 end 142 T([pos_k; pos_kminus]) = []; 143 rho([pos_k; pos_kminus]) = []; 144 145 % Casos en que sea necesario eliminar una arista e introducir otra (si esta no 146 % va en la dirección del segmento ij_{j0}) 147 elseif ~isnan(pos_k) && isnan(pos_kminus) 148 149 % Vector director del segmento ij_{j0} 150 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 151 152 % Vector t 153 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 154 155 if (abs(t_0(1) - v_0(1)) < 1e-12 && abs(t_0(2) - v_0(2)) < 1e-12) || ... 156 (abs(t_0(1) + v_0(1)) < 1e-12 && abs(t_0(2) + v_0(2)) < 1e-12) % Se ha 157 % aplicado una cierta tolerancia 158 % Caso en que solo se suprime una arista 159 Tinv(k) = NaN; Tinv(T(pos_k + 1:length(T))) = Tinv(T(pos_k + 1: ... 160 length(T))) - 1; 161 T(pos_k) = []; 162 rho(pos_k) = []; 163 else 164 % Caso en que se sustituye una arista por la otra 165 Tinv(k) = NaN; Tinv(kminus) = pos_k; 166 T(pos_k) = kminus; 167 end 168 169 % Caso análogo al anterior 170 elseif isnan(pos_k) && ~isnan(pos_kminus) 171 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 172 u = V_Mercator(k+1,:) - V_Mercator(k,:); u_0 = u/norm(u); 173 174 if (abs(u_0(1) - v_0(1)) < 1e-12 && abs(u_0(2) - v_0(2)) < 1e-12) || ...
D.3 Algoritmo de Lee 97 175 (abs(u_0(1) + v_0(1)) < 1e-12 && abs(u_0(2) + v_0(2)) < 1e-12) 176 Tinv(kminus) = NaN; Tinv(T(pos_kminus + 1:length(T))) = ... 177 Tinv(T(pos_kminus + 1:length(T))) - 1; 178 T(pos_kminus) = []; 179 rho(pos_kminus) = []; 180 else 181 Tinv(kminus) = NaN; Tinv(k) = pos_kminus; 182 T(pos_kminus) = k; 183 end 184 185 % Caso en que hay que añadir dos aristas (si la dirección de estas no va en la 186 % del segmento ij_{j0}) 187 else 188 189 % Vector director del segmento ij_{j0} 190 v = (WPs_Mercator(j(j0),:) - WPs_Mercator(i,:))’; v_0 = v/norm(v); 191 192 % Vectores u y t 193 u = V_Mercator(k + 1,:) - V_Mercator(k,:); u_0 = u/norm(u); 194 t = V_Mercator(kminus,:) - V_Mercator(k,:); t_0 = t/norm(t); 195 196 % Distancia al nodo i 197 rho_k0 = norm(v); 198 199 % Definición de las aristas a introducir en T 200 if (abs(u_0(1) - v_0(1)) < 1e-12 && abs(u_0(2) - v_0(2)) < 1e-12) || ... 201 (abs(u_0(1) + v_0(1)) < 1e-12 && abs(u_0(2) + v_0(2)) < 1e-12) 202 newedges = kminus; % Solo se añadirá una arista 203 204 elseif (abs(t_0(1) - v_0(1)) < 1e-12 && abs(t_0(2) - v_0(2)) < 1e-12) ... 205 || (abs(t_0(1) + v_0(1)) < 1e-12 && abs(t_0(2) + v_0(2)) < 1e-12) 206 newedges = k; % Solo se añadirá una arista 207 208 else 209 210 %Se añadirán dos aristas 211 rho_k0 = [rho_k0; rho_k0]; 212 213 % Las dos aristas se ordenan a través del producto escalar 214 udotv = dot(u_0,v); 215 tdotv = dot(t_0,v); 216 217 if udotv > tdotv 218 newedges = [kminus; k]; 219 else 220 newedges = [k; kminus]; 221 end 222 223 end 224 225 % Introducción en T y actualización de rho y Tinv 226 if isempty(T) 227 T = newedges; 228 rho = rho_k0; 229 Tinv(T(1:length(newedges))) = 1:length(newedges); 230 else 231 [pos_newedges, rho, T] = BinarySearchAndInsertion(rho, rho_k0, T, ... 232 newedges, V_Mercator, WPs_Mercator(i,:), WPs_Mercator(j(j0),:)); 233 Tinv(T(pos_newedges)) = pos_newedges; 234 Tinv(T(pos_newedges(length(pos_newedges)) + 1:length(T))) = ... 235 Tinv(T(pos_newedges(length(pos_newedges)) + 1:length(T))) + ... 236 length(pos_newedges);