ADVERTIMENT. L'accés als continguts d'aquesta tesi doctoral i la seva utilització ha de respectar els drets de la persona autora. Pot ser utilitzada per a consulta o estudi personal, així com en activitats o materials d'investigació i docència en els termes establerts a l'art. 32 del Text Refós de la Llei de Propietat Intel·lectual (RDL 1/1996). Per altres utilitzacions es requereix l'autorització prèvia i expressa de la persona autora. En qualsevol cas, en la utilització dels seus continguts caldrà indicar de forma clara el nom i cognoms de la persona autora i el títol de la tesi doctoral. No s'autoritza la seva reproducció o altres formes d'explotació efectuades amb finalitats de lucre ni la seva comunicació pública des d'un lloc aliè al servei TDX. Tampoc s'autoritza la presentació del seu contingut en una finestra o marc aliè a TDX (framing). Aquesta reserva de drets afecta tant als continguts de la tesi com als seus resums i índexs.

ADVERTENCIA. El acceso a los contenidos de esta tesis doctoral y su utilización debe respetar los derechos de la persona autora. Puede ser utilizada para consulta o estudio personal, así como en actividades o materiales de investigación y docencia en los términos establecidos en el art. 32 del Texto Refundido de la Ley de Propiedad Intelectual (RDL 1/1996). Para otros usos se requiere la autorización previa y expresa de la persona autora. En cualquier caso, en la utilización de sus contenidos se deberá indicar de forma clara el nombre y apellidos de la persona autora y el título de la tesis doctoral. No se autoriza su reproducción u otras formas de explotación efectuadas con fines lucrativos ni su comunicación pública desde un sitio ajeno al servicio TDR. Tampoco se autoriza la presentación de su contenido en una ventana o marco ajeno a TDR (framing). Esta reserva de derechos afecta tanto al contenido de la tesis como a sus resúmenes e índices.

WARNING. Access to the contents of this doctoral thesis and its use must respect the rights of the author. It can be used for reference or private study, as well as research and learning activities or materials in the terms established by the 32nd article of the Spanish Consolidated Copyright Act (RDL 1/1996). Express and previous authorization of the author is required for any other uses. In any case, when using its content, full name of the author and title of the thesis must be clearly indicated. Reproduction or other forms of for profit use or public communication from outside TDX service is not allowed. Presentation of its content in a window or frame external to TDX (framing) is not authorized either. These rights affect both the content of the thesis and its abstracts and indexes.

´ ˜ UNIVERSITAD POLITECNICA DE CATALUNA

Programa de Doctorado: ´ AVANZADA Y ROBOTICA ´ AUTOMATIZACION

Tesis Doctoral

Planificaci´ on de Movimientos en Entornos Din´ amicos o Inciertos Mediante la Coordinaci´ on de M´ etodos Aleatorios de B´ usqueda y Funciones Arm´ onicas Pedro I˜ niguez Galbete

Director: Jan Rosell Gratac`os

Instituto de Organitzaci´on y Control de Sistemas Industriales Febrero del 2012

Resumen En los m´etodos planificadores de trayectorias basados en funciones potenciales, la utilizaci´on de las funciones arm´onicas tiene la importante propiedad de no presentar m´ınimos locales. Sin embargo, la creaci´on de planificadores basados en estas funciones arm´onicas se ha encontrado con serias dificultades, sobre todo cuando el n´ umero de grados de libertad es elevado. Por este motivo, esta tesis realiza inicialmente un estudio de las propiedades m´as relevantes de dichas funciones arm´onicas; destacando aquellas que han sido la causa de su reducida aplicaci´on en la generaci´on de trayectorias. Al mismo tiempo, el resultado de este estudio sirve de base para la proposici´on de m´etodos compensatorios que permitan reducir las propiedades negativas de las funciones arm´onicas, como funciones potenciales aplicables a la generaci´on de movimientos en rob´otica. Despu´es se considera los m´etodos num´ericos de c´alculo de las funciones arm´onicas, as´ı como el coste computacional de los mismos. Con el objetivo de reducir el tiempo de c´alculo, esta tesis propone una discretizaci´on jer´arquica y un m´etodo eficiente de etiquetado de celdas. Por su parte, dicha discretizaci´on jer´arquica, se va realizando progresivamente mediante muestreo aleatorio y descomposici´on de celdas, lo que genera un escenario parcialmente conocido que, sin embargo, permitir´a en cierto n´ umero de casos encontrar la trayectoria buscada. Por lo tanto, esta propuesta reduce dr´asticamente el n´ umero de puntos de c´alculo y, por consiguiente, el tiempo de computaci´on. La tesis completa la propuesta de un planificador combinando las t´ecnicas de muestreo con el c´alculo de funciones arm´onicas mediante un m´etodo de exploraci´on aleatorio conducido (PHM), aplicado a un espacio de configuraciones discretizado jer´arquicamente sobre el que se va recalculando la funci´on arm´onica. De esta forma la exploraci´on se gu´ıa hacia zonas m´as prometedoras, intentando obtener la soluci´on por fases.

´Indice general 1. Introducci´ on 1.1. Marco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1. Espacio de Configuraciones . . . . . . . . . . . . . . 1.2.1.1. C´alculo del espacio de configuraciones . . 1.2.1.2. Discretizaci´on y descomposici´on en celdas 1.2.1.3. Muestreo del espacio de configuraciones . 1.2.2. Planificaci´on de movimientos . . . . . . . . . . . . 1.2.2.1. M´etodos completos . . . . . . . . . . . . . 1.2.2.2. M´etodos basados en muestreo . . . . . . . 1.2.2.3. M´etodos para entornos din´amicos . . . . . 1.3. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

1 2 3 3 3 5 6 7 7 9 11 11

2. Planificaci´ on con funciones potenciales arm´ onicas, subarm´ onicas y superarm´ onicas 15 2.1. Introducci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2. Propiedades afines a la generaci´on de funciones potenciales para la planificaci´on de movimientos en rob´otica . . . . . . . . 18 2.2.1. Ausencia de m´ınimos locales . . . . . . . . . . . . . . . 19 2.2.1.1. An´alisis en dos dimensiones . . . . . . . . . . 19 2.2.1.2. An´alisis en n dimensiones . . . . . . . . . . . 22 2.2.2. Concavidad, convexidad de las funciones subarm´onicas y superarm´onicas . . . . . . . . . . . . . . . . . . . . . 24 2.2.3. L´ıneas de corriente y superficies equipotenciales . . . . 25 2.2.4. Superposici´on lineal . . . . . . . . . . . . . . . . . . . . 26 2.2.5. Curvatura de las funciones arm´onicas . . . . . . . . . . 26 2.2.6. Estabilidad de la soluci´on . . . . . . . . . . . . . . . . 27 2.3. Problem´atica inherente al m´etodo . . . . . . . . . . . . . . . . 27 2.3.1. Efecto embudo y funciones arm´onicas . . . . . . . . . . 27 2.3.2. Condiciones de contorno, condiciones de solubilidad y balance global del flujo del gradiente . . . . . . . . . . 29 2.4. C´alculo anal´ıtico y obtenci´on num´erica de soluciones . . . . . 31 2.4.1. Modelo base. An´alisis de la influencia de las condiciones de contorno mediante las Funciones de Green . . . 32 2.4.2. Obtenci´on de soluciones mediante c´alculo num´erico . . 33 2.5. Aportaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 i

ii

´INDICE

3. Alternativas de moldeado de las funciones arm´ onicas para su uso en la planificaci´ on de movimientos 37 3.1. Modelo gen´erico . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2. Condiciones de contorno de Dirichlet. Obtenci´on de la funci´on de Green mediante el m´etodo de las cargas im´agenes, aplicado a funciones de dos variables . . . . . . . . . . . . . . . . . . . 39 3.2.1. Funci´on de Green para todo el espacio . . . . . . . . . 39 3.2.2. Ejemplos basados en condiciones de contorno homog´eneas: Rect´angulo y circunferencias obst´aculo . . . . . . . . . 40 3.2.3. Contorno infinito, funci´on delta y obst´aculo circular . . 46 3.3. Moldeado de la funci´on potencial mediante la aplicaci´on de condiciones de contorno no uniformes de Dirichlet . . . . . . . 49 3.3.1. Dominio en forma de franja . . . . . . . . . . . . . . . 49 3.3.2. Condiciones de contorno proporcional con la distancia sobre el contorno de un obst´aculo circular . . . . . . . 52 3.4. Moldeado de la funci´on potencial mediante el uso de condiciones de contorno mixtas . . . . . . . . . . . . . . . . . . . . . . 52 3.4.1. Obtenci´on de la funci´on de Green utilizando las series de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.2. Propuesta de moldeado . . . . . . . . . . . . . . . . . . 54 3.5. Moldeado de la funci´on potencial mediante la utilizaci´on de las funciones superarm´onicas y subarm´onicas . . . . . . . . . . 54 3.5.1. Contorno infinito, funci´on delta y obst´aculo circular . 54 3.5.2. Contorno infinito, funci´on delta, obst´aculos y parche en torno del punto objetivo. . . . . . . . . . . . . . . . 57 3.5.3. Modelizaci´on de los obst´aculos mediante coronas circulares. . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.6. Aportaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4. Planificaci´ on de movimientos en entornos din´ amicos y cambiantes 4.1. Estudio anal´ıtico para condiciones iniciales y funci´on obst´aculo 4.2. Control de la configuraci´on y de la trayectoria de un robot en un entorno con obst´aculos m´oviles. . . . . . . . . . . . . . . . 4.3. Distribuci´on espacial del cambio de la funci´on potencial debido a un cambio local . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Aportaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. C´ alculo num´ erico de las funciones arm´ onicas, superarm´ onicas y subarm´ onicas 5.1. Discretizaci´on uniforme . . . . . . . . . . . . . . . . . . . . . . 5.1.1. Ecuaci´on en diferencias . . . . . . . . . . . . . . . . . . 5.1.2. Condiciones de contorno de Dirichlet y de Neuman . . 5.1.3. M´etodos iterativos . . . . . . . . . . . . . . . . . . . . 5.1.3.1. M´etodo iterativo de Jacobi . . . . . . . . . . 5.1.3.2. M´etodo iterativo de Gaus-Seidel . . . . . . . .

65 65 69 71 77

79 80 80 82 83 84 85

´INDICE

5.2. 5.3.

5.4. 5.5.

iii 5.1.3.3. M´etodo iterativo S.O.R. (Sucesive Over Relaxation) . . . . . . . . . . . . . . . . . . . 5.1.4. C´alculo en multirresoluci´on . . . . . . . . . . . . . . Discretizaci´on jer´arquica adaptativa . . . . . . . . . . . . . . 5.2.1. Esquema num´erico de identificaci´on de celdas . . . . Multirresoluci´on Jer´arquica Adaptativa . . . . . . . . . . . . 5.3.1. Interpolaci´on . . . . . . . . . . . . . . . . . . . . . . 5.3.2. Mapeado . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.3. Algoritmos . . . . . . . . . . . . . . . . . . . . . . . . 5.3.3.1. Algoritmo de c´alculo de la funci´on potencial arm´onica del C-espacio . . . . . . . . . . . . 5.3.3.2. Algoritmo de obtenci´on de la trayectoria . . Obtenci´on de resultados . . . . . . . . . . . . . . . . . . . . Aportaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

86 88 90 91 93 95 96 97

. . . .

97 98 100 104

6. Uso de las funciones arm´ onicas en espacios de configuraciones explorados mediante t´ ecnicas de muestreo 105 6.1. Planificador PHM b´asico . . . . . . . . . . . . . . . . . . . . . 106 6.1.1. Selecci´on aleatoria de una celda ponderada por su tama˜ no: sorteo para la exploraci´on de celdas . . . . . . . 107 6.1.2. Verificaci´on y clasificaci´on de celdas . . . . . . . . . . . 109 6.1.3. Algoritmo del planificador de trayectorias. Exploraci´on pasiva de trayectorias . . . . . . . . . . . . . . . . . . . 110 6.1.4. Estudio de la completitud del m´etodo propuesto . . . . 112 6.1.4.1. Modelo de poblaci´on de celdas en G . . . . . 113 6.1.4.2. Probabilidad de selecci´on condicionada por el tama˜ no de las celdas . . . . . . . . . . . . . . 115 6.1.4.3. Probabilidad de fallo . . . . . . . . . . . . . . 116 6.1.4.4. Validaci´on . . . . . . . . . . . . . . . . . . . . 117 6.2. Alternativas de mejora del m´etodo b´asico . . . . . . . . . . . . 118 6.2.1. Exclusi´on de celdas del sorteo de exploraci´on . . . . . . 118 6.2.2. Exploraci´on pasiva de trayectorias: Puntos de Apoyo y Funci´on Potencial Auxiliar . . . . . . . . . . . . . . . . 119 6.2.3. Resoluci´on m´axima ajustable . . . . . . . . . . . . . . 120 6.2.4. Aplicaci´on de condiciones de contorno de Neuman en los obst´aculos . . . . . . . . . . . . . . . . . . . . . . . 120 6.3. Algoritmo del planificador PHM avanzado . . . . . . . . . . . 121 6.4. Exploraci´on basada en la detecci´on de colisi´on: Planificador PHM color . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.4.1. Secuencia determinista de muestreo . . . . . . . . . . . 125 6.4.2. Clasificaci´on probabil´ıstica de celdas . . . . . . . . . . 125 6.4.3. Muestreo de celdas . . . . . . . . . . . . . . . . . . . . 127 6.4.4. Algoritmo del planificador PHM color . . . . . . . . . . 129 6.5. Reducci´on del tiempo de computaci´on respecto de una discretizaci´on uniforme . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.6. Aportaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

iv

´INDICE

7. Conclusiones 133 7.1. Publicaciones realizadas con el desarrollo de la tesis . . . . . . 134 7.2. Trabajos futuros . . . . . . . . . . . . . . . . . . . . . . . . . 135 A. Demostraciones 137 A.1. Principio de localidad de una funci´on . . . . . . . . . . . . . . 137 ´ A.2. Area e integraci´on de una hiperesfera . . . . . . . . . . . . . . 140 A.3. Propiedades de las funciones arm´onicas, superarm´onicas y subarm´onicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 A.3.1. Concavidad, convexidad de las funciones arm´onicas y superarm´onicas . . . . . . . . . . . . . . . . . . . . . . 143 A.3.2. Curvatura y optimizaci´on de las trayectorias obtenidas con funciones arm´onicas . . . . . . . . . . . . . . . . . 146 A.3.3. Condiciones de contorno de Dirichlet . . . . . . . . . . 148 A.3.4. Condiciones de contorno de Neuman . . . . . . . . . . 150 A.3.5. Condiciones de contorno mixtas (Robin o de radiaci´on) 150 A.3.6. Estabilidad de la soluci´on . . . . . . . . . . . . . . . . 151 A.4. Teorema de Green . . . . . . . . . . . . . . . . . . . . . . . . 152 A.5. M´etodos iterativos . . . . . . . . . . . . . . . . . . . . . . . . 153 A.5.1. M´etodo iterativo de Jacobi . . . . . . . . . . . . . . . . 153 A.5.2. M´etodo iterativo de Gaus-Seidel . . . . . . . . . . . . . 157 A.5.3. M´etodo iterativo S.O.R . . . . . . . . . . . . . . . . . . 159 A.5.4. M´etodo iterativo ponderado de Jacobi . . . . . . . . . 161

Cap´ıtulo 1 Introducci´ on El t´ermino rob´otica introducido por Isaac Asimov corresponde a un concepto idealizado que supera las pragm´aticas y estandarizadas definiciones, en las cuales un robot industrial es un manipulador que funciona m´as como una m´aquina repetitiva que como un sistema vers´atil e inteligente. El desarrollo tecnol´ogico de la micro electr´onica y la micro inform´atica junto con el avance en disciplinas tales como la visi´on artificial, elementos h´apticos, integraci´on sensorial e inteligencia artificial han permitido progresar y crear nuevas aplicaciones en rob´otica m´edica, rob´otica m´ovil, robots humanoides y microrob´otica. Por otro lado, la necesidad de disponer de robots industriales m´as vers´atiles y aut´onomos que permitan operaciones de ensamblado complejas, ha impulsado la investigaci´on y desarrollo de subsistemas inteligentes generadores de secuencias y de las trayectorias de configuraciones correspondientes: Se precisa programar a nivel de tarea y no a nivel de movimientos elementales, ahorrar tiempo de programaci´on y realizar aplicaciones para entornos cambiantes e imprecisos. Al mismo tiempo, las nuevas metodolog´ıas y t´ecnicas desarrolladas permiten un nuevo enfoque en otras disciplinas, ofreciendo nuevas posibilidades en el an´alisis y dise˜ no, y extendiendo el concepto de rob´otica hacia una idea m´as global. Gen´ericamente, un robot es un cuerpo m´ovil que interacciona con otros dentro de un espacio f´ısico. Puede estar formado por uno o varios elementos, ligados o fijos, r´ıgidos o flexibles; pudiendo concretarse desde una pieza para ensamblar hasta un robot manipulador industrial. En un sentido amplio, un sistema robotizado aut´onomo, compuesto por uno o varios brazos manipuladores, ha de ser capaz de realizar diversas operaciones previamente definidas, para lo cual necesita un sistema inteligente constituido por varios subsistemas interrelacionados jer´arquicamente (figura 1.1): el Planificador de tareas, el Planificador de movimientos y el Generador de trayectorias. Este trabajo se centra en el Planificador de movimientos y desarrolla un m´etodo nuevo denominado Probabilistic Harmonic function based Method; el cual mejora la eficiencia temporal de computaci´on respecto de los m´etodos de campo de potencial que utilizan las funciones arm´onicas, presentando un 1

´ CAP. 1. INTRODUCCION

2

Figura 1.1: Jerarqu´ıa de funciones en un sistema robotizado.

comportamiento completo; probabil´ısticamente y en resoluci´on.

1.1.

Marco

En un Sistema robotizado aut´onomo, el subsistema Planificador de trayectorias est´a supeditado, jer´arquicamente, al subsistema Planificador de tareas; el cual convierte un objetivo global en operaciones elementales de desplazamiento, ensamblado, mecanizado y otros. El subsistema Planificador de movimientos debe ser capaz de generar la trayectoria entre una configuraci´on inicial y otra final, autom´aticamente, sin colisiones con los objetos del entorno pero con posibles desplazamientos correctivos y contactos adaptativos de ensamblado. Esta funci´on debe de realizarse con eficiencia temporal y de forma completa; es decir, si existe una soluci´on debe de encontrarla en un tiempo m´ınimo y si no, debe de indicar la imposibilidad del movimiento. Por su parte, el Generador de trayectorias determina la velocidades de los movimientos para cada una de las variables o articulaciones. Los movimientos desarrollados por el robot en la ejecuci´on de una tarea se clasifican en: movimientos gruesos de acercamiento o aproximaci´on, en los cuales no est´a permitido ning´ un contacto con los obst´aculos, (gross-motion planning) y movimientos finos (fine-motion planning) donde los contactos de acomodaci´on mediante el control de fuerzas, en la consecuci´on de operaciones de ensamblado, son parte esencial y donde la incertidumbre constituye un problema central. Generalmente la planificaci´on de movimientos se realiza en el espacio de configuraciones, donde la parametrizaci´on del robot permite representar una configuraci´on del mismo en un punto, pudiendo ser libre o de colisi´on; de forma que todos los puntos libres definen el subespacio libre y todos los puntos de colisi´on el subespacio obst´aculo. Dichos subespacios constituyen el entorno sobre el cual se planificar´an los movimientos. En algunos casos sencillos de robots m´oviles desplaz´andose en el plano y algunos manipuladores, se realiza la planificaci´on de movimientos directamente sobre el espacio f´ısico.

´ CAP. 1. INTRODUCCION

1.2.

3

Antecedentes

El estudio del estado del arte se realiza considerando que la planificaci´on de movimientos en rob´otica puede separarse en dos etapas [68]: 1. Determinar el espacio de configuraciones del robot. 2. Encontrar la trayectoria libre de colisiones entre un punto inicial y otro final.

1.2.1.

Espacio de Configuraciones

En rob´otica la planificaci´on de movimientos se realiza usualmente en el espacio de configuraciones [68] (C-space), donde la configuraci´on del robot se transforma en un punto mediante la parametrizaci´on de la posici´on o configuraci´on f´ısica del robot, esta importante propiedad reduce un movimiento en el espacio f´ısico a una l´ınea en el espacio de configuraciones. El n´ umero de dimensiones del espacio de configuraciones coincide con el n´ umero de grados de libertad del robot. As´ı, un cuerpo r´ıgido en el espacio 3D tiene seis grados de libertad (tres de posici´on y tres de orientaci´on), por lo que el espacio de configuraciones correspondiente ser´a de seis dimensiones y si se disponen de n cuerpos ser´a de 6n. Tambi´en, el espacio de configuraciones correspondiente a un manipulador, formado por una cadena de n cuerpos r´ıgidos ligados por n articulaciones, dispone de n dimensiones y una determinada configuraci´on del mismo se transforma en un punto definido por n coordenadas [62]. En general, en el espacio f´ısico de trabajo tendremos un n´ umero determinado de cuerpos cuyas posiciones pueden ser libres, en contacto o prohibidas, de manera que se transforman en el espacio de configuraciones en C-obst´aculos cuyos contornos definen las superficies de contactos, la zona interior ser´a prohibida y la exterior libre. La determinaci´on del espacio de configuraciones no es trivial [22] y ha sido motivo de investigaci´on durante las u ´ltimas d´ecadas. La dificultad de transformaci´on depende del n´ umero de dimensiones, de si los objetos son convexos o c´oncavos y del n´ umero de v´ertices de los mismos, asumiendo que son poli´edricos. Previamente a la obtenci´on del espacio de configuraciones es necesario modelar los objetos existentes en el espacio de trabajo, Wise y Bowyer [94] clasifican las posibles estrategias en la obtenci´on del espacio de configuraciones en funci´on de la precisi´on en el modelado de los objetos y en la precisi´on de desarrollo del espacio de configuraciones. 1.2.1.1.

C´ alculo del espacio de configuraciones

Se han presentado varios m´etodos para computar el espacio de configuraciones [94]: Computaci´ on anal´ıtica de las Superficies de Contacto: La obtenci´on de expresiones que definan anal´ıticamente las superficies de contacto es especialmente necesaria en movimientos de precisi´on; donde es preciso controlar movimientos de acomodaci´on en contacto para la consecuci´on de operaciones de ensamblado. Para ello es fundamental la clasificaci´on de contactos

4

´ CAP. 1. INTRODUCCION

b´asicos realizada por Lozano-P´erez [68], las condiciones algebraicas de los mismos, en un determinado rango, definen las superficies de contacto simple. La intersecci´on de estas superficies determina las aristas de contacto doble y la intersecci´on de tres aristas un punto de contacto triple. Lozano-P´erez en [70] realiza una discretizaci´on de las variables generalizadas, con proyecciones recursivas hasta reducir el problema a dos dimensiones donde un objeto libre se desliza en contacto sobre el contorno poligonal del obst´aculo. Para este caso Brost [17] describe un algoritmo que computa el espacio de configuraciones con una descripci´on topol´ogica de las superficies de los C-obst´aculos. Rosell, Basa˜ nez y Su´arez [86] construyen el espacio de configuraciones para el caso de un pol´ıgono m´ovil 2D con desplazamiento y rotaci´on, en contacto con otro pol´ıgono fijo. Las superficies 3D desarrolladas en el espacio de configuraciones se proyectan sobre el plano desplazamiento, reduciendo una dimensi´on y facilitando su utilizaci´on en la planificaci´on de movimientos para ensamblaje en 2D. Hwang [46] describe un m´etodo para obtener las ecuaciones de los bordes de los obst´aculos en el espacio de configuraciones basado en la aproximaci´on de los obst´aculos como poliedros de caras triangulares y los elementos del brazo manipulador como segmentos, la intersecci´on de estos dos elementos define una ecuaci´on de la que se obtiene las expresiones anal´ıticas correspondientes. Maciejewski [71] realiza la computaci´on anal´ıtica solo de los bordes de aquellos obst´aculos cuya colisi´on es posible y determina las tangentes de las curvas, lo cual ayuda al conocimiento topol´ogico del espacio de configuraciones. La complejidad depende del n´ umero de v´ertices de obst´aculos y del n´ umero de grados de libertad, si se trata de un robot s´olido que se desplaza y gira o se trata de un brazo manipulador en 3D. Construcci´ on del mapa de bits del espacio de configuraciones: Los obst´aculos y espacio libre son representados en un mapa de bits por unos y ceros respectivamente, lo cual reduce dr´asticamente el tiempo de comprobaci´on de colisi´on. Este mapa de bits puede ser obtenido discretizando las expresiones anal´ıticas de los procedimientos anteriores, o bien observando que, cuando el robot es un objeto r´ıgido con movimiento de translaci´on, el espacio de configuraciones es la convoluci´on entre el espacio de trabajo y los obst´aculos. Kavraki propone en [50] el uso de la Transformada R´apida de Fourier para computar esta convoluci´on en el caso de un robot m´ovil y Curto y Moreno [29] generalizan este m´etodo para el caso de robots manipuladores. Lozano-P´erez [68] muestra para el caso de un robot m´ovil con desplazamiento sin rotaci´on, la obtenci´on de los C-obst´aculos por el crecimiento del perfil de los obst´aculos a˜ nadiendo el del robot en contacto con el contorno, usando la suma Minkowski. Newman y Branicky [82] identifican formas elementales de objetos, en el espacio de trabajo, que son f´acilmente transformables al espacio de configuraciones. Estas transformaciones de formas elementales son memorizadas como mapa de bits y combinadas para la obtenci´on de formas m´as complejas. El tiempo de computaci´on del mapa de bits, as´ı como la memoria requerida, son funci´on de su tama˜ no y resoluci´on. As´ı, en canales estrechos y formas complicadas se requiere una resoluci´on elevada y, si el n´ umero de grados de libertad es tambi´en elevado, el tiempo de computaci´on y la memoria de almacenamiento de variables crecer´an exponencialmente.

´ CAP. 1. INTRODUCCION 1.2.1.2.

5

Discretizaci´ on y descomposici´ on en celdas

Son dos m´etodos conceptualmente distintos utilizados con fines distintos: mientras que una discretizaci´on es una cuantificaci´on realizada con la finalidad de realizar c´alculos num´ericos de funciones, la descomposici´on en celdas es una partici´on conexa de regiones a las que se le hace corresponder un grafo. La discretizaci´on puede ser regular o no regular. Si queremos un conocimiento exhaustivo del espacio de configuraciones deberemos de realizar una discretizaci´on del mismo con un paso igual o menor a la precisi´on requerida. Utilizando el m´etodo 2n -tree ([102], [23], [43] y [47]) obtendremos una discretizaci´on no uniforme y jerarquizada, que se adapta a los contornos de los obst´aculos con las celdillas de mayor resoluci´on; utiliza celdas gruesas en el espacio libre y s´olido de los obst´aculos, reduciendo el n´ umero de puntos para el c´alculo. Esta discretizaci´on jerarquizada genera celdas contenidas en el espacio libre (blancas), celdas contenidas en el interior del obst´aculo (negras) y celdas que contienen el borde del obst´aculo (grises), las cuales son susceptibles de ser nuevamente particionadas por el procedimiento 2n -tree Blanco et al [9] presentan un algoritmo que eval´ ua el espacio de configuraciones de un robot m´ovil con dos grados de libertad. Se basa en la convoluci´on discreta del espacio de trabajo y el robot, la cual se realiza sobre una partici´on jer´arquica quad-tree, empezando por el nivel m´as bajo hasta llegar a la m´axima resoluci´on. En cada paso se determina una ventana que contiene obst´aculo, lo que reduce el n´ umero de puntos para el c´alculo. La reducci´on de memoria y tiempo de c´omputo es mayor cuanto menor es el n´ umero de obst´aculos. La descomposici´on en celdas puede ser exacta o aproximada seg´ un su uni´on coincida exactamente con el espacio libre o est´e contenida en ´el. La descomposici´on exacta se realiza mediante celdas no solapadas y definidas por puntos significativos de los C-obst´aculos. Si el espacio de configuraciones es de dos dimensiones y los C-obst´aculos son pol´ıgonos, las celdas ser´an tri´angulos, rect´angulos o trapecios con sus v´ertices coincidentes con los v´ertices de los C-obst´aculos. Similarmente se procede para el caso de un espacio de configuraciones 3D, apareciendo en este caso celdas de formas m´as complejas. La descomposici´on aproximada deja una cierta cantidad de espacio libre sin adjudicaci´on de celda, por lo que los m´etodos planificadores basados en esta descomposici´on no son completos; es decir, es posible que no encuentren una trayectoria a pesar de que esta exista. Por este motivo muchos m´etodos de descomposici´on aproximada utilizan celdas de tama˜ no sucesivamente adaptable a los contornos de los C-obst´aculos. Kambhampati y Davis [48] utilizan una descomposici´on en multiresoluci´on. Tambi´en en la descomposici´on aproximada el m´etodo es ampliamente utilizado. Obtendremos una partici´on en celdas de resoluci´on adaptable a los contornos de los obst´aculos. Esta partici´on jerarquizada permite un etiquetado de celdas relacionadas por niveles; lo que facilita la formaci´on del ´arbol de b´ usqueda correspondiente. As´ı, Harabor y Botea en [42] realizan una partici´on del espacio libre en celdas de diferentes tama˜ nos que, interconectadas

6

´ CAP. 1. INTRODUCCION

transversalmente, constituyen un grafo de agentes de tama˜ nos heterog´eneos. Un algoritmo de busqueda sobre dicho grafo obtendr´a la trayectoria deseada.Tambi´en Zhu y Latombe en [103] consideran la descomposici´on jer´arquica, de celdas aproximadas, del espacio de configuraciones. Con dichas celdas obtienen el grafo de conectividad del espacio libre. Entonces dise˜ nan nuevo algoritmo de b´ usqueda jer´arquica, con un mecanismo para el registro de las condiciones de fallo. Por su parte, Cai y Ferrari en [19] desarrollan un m´etodo de descomposici´on aproximada de celdas; en el que los obst´aculos, los objetivos y la plataforma de sensores est´an representados como subconjuntos cerrados y acotados, en un espacio de trabajo eucl´ıdeo. El m´etodo construye un grafo de conectividad con las celdas, que se poda y se transforma en un ´arbol en el que se puede calcular una estrategia ´optima de b´ usqueda. • En los m´etodos jer´arquicos de partici´on son fundamentales aquellas funciones detectoras de colisi´on que determinan la distancia m´ınima a los obst´aculos [65], [92]. En algunos casos de particiones jer´arquicas, sobre las que se calcula una funci´on potencial, es importante poder disponer de una funci´on detectora de colisi´on que nos reporte tambi´en la distancia de penetraci´on para los casos de colisi´on. Tanto la distancia m´ınima al obst´aculo como la distancia de penetraci´on se computan en el espacio f´ısico, por lo que deben de transformarse al espacio de configuraciones. 1.2.1.3.

Muestreo del espacio de configuraciones

La eficiencia computacional de las funciones detectoras de colisi´on dan pie al enfoque basado en muestreo. Siguiendo una estrategia: aleatoria, determinista o combinaci´on de ambas se exploran una serie de puntos del espacio de configuraciones y con la informaci´on obtenida, de colisi´on o libertad, de estos puntos es posible plantear m´etodos algor´ıtmicos que deduzcan trayectorias libres de colisi´on. Cuando el n´ umero de grados de libertad es elevado el tiempo de computaci´on empleado en el c´alculo de las superficies de contacto, la obtenci´on del mapa de bits del espacio de configuraciones o la discretizaci´on del mismo, es tambi´en muy elevado. Por otro lado, no es necesario tener un conocimiento completo del espacio de configuraciones para obtener una trayectoria de movimiento, libre de colisiones, entre dos puntos. Simplemente con el conocimiento de un determinado n´ umero de puntos libres y de obst´aculo es posible obtener, en determinadas condiciones, las trayectorias buscadas. Con la finalidad de buscar posiciones libres y su conectividad, las funciones detectoras de colisi´on debe determinar si un punto en el espacio de configuraciones corresponde al espacio obst´aculo o al espacio libre. Se utilizan en espacios de configuraciones con un n´ umero elevado de grados de libertad, y como funci´on recurrente de planificadores de movimientos basados en m´etodos probabil´ısticos. En estos casos, el tiempo de computaci´on de la trayectoria depende del tiempo de computaci´on de la funci´on y del n´ umero de veces que sea invocada. Con este objetivo, Mirtich [79] describe y analiza varios algoritmos de de-

´ CAP. 1. INTRODUCCION

7

tecci´on utilizados en la planificaci´on de movimientos, tanto en movimientos gruesos como finos. Destaca conceptos generales como coherencia y localidad, as´ı como la generalizada estrategia de detecci´on de colisi´on en dos fases: fase ancha y fase estrecha. Mientras que Lin y Canny [65] desarrollan un algoritmo detector de colisiones muy eficiente para el caso de c´alculos repetitivos, con objetos en movimiento. Realizan una descomposici´on de los objetos en formas b´asicas, sobre las que definen una serie de situaciones sencillas: par de v´ertices, v´ertice y arista, v´ertice y cara, par de aristas, arista y cara, dos caras. Sucesivamente va comprobando cada par de situaciones entre dos objetos, hasta encontrar la distancia m´ınima entre ellos. Mart´ınez-Salvador y del Pobil [72] realizan una descomposici´on de los objetos en 3D mediante esferas. La simplicidad aportada por esta partici´on reduce el tiempo de computaci´on en la detecci´on de colisiones (P´erez y del Pobil [83]) y en el seguimiento de objetos en movimiento. Por su parte, Kim, Lin y Manocha [57] presentan un algoritmo incremental para estimar la distancia de penetraci´on entre objetos convexos en 3D. Se basa en el c´alculo de las sumas de Minkowski y requiere para su buen funcionamiento alta coherencia espacial y de movimiento. El conocimiento de la distancia a los obst´aculos es tambi´en utilizado en el m´etodo propuesto por Gilbert y Johnson en [36], donde la idea principal consiste en expresar la fuerza de repulsi´on de los obst´aculos en funci´on de la distancia entre las partes que pueden colisionar.

1.2.2.

Planificaci´ on de movimientos

Los m´etodos de planificaci´on de movimientos existentes se pueden clasificar en funci´on de su capacidad para encontrar una soluci´on, lo que a su vez est´a relacionado con el tiempo empleado para su c´alculo. As´ı, existen dos m´etodos principales que se describen a continuaci´on: M´etodos completos y M´etodos basados en muestreo. Dentro de los m´etodos basados en muestreo hay que distinguir, a su vez, aquellos que utilizan una ley determinista para generar las configuraciones de muestreo, y los que utilizan procedimientos aleatorios. Adem´as, hay que distinguir los m´etodos que consideran entornos est´aticos de aquellos que consideran entornos din´amicos, con la posibilidad de uno o varios obst´aculos en movimiento. 1.2.2.1.

M´ etodos completos

Los m´etodos completos son aquellos que encuentran una soluci´on si esta existe, y en el caso posible de que no exista lo indica mediante la devoluci´on de una variable. Si bien esta opci´on es deseable el coste computacional correspondiente es elevado, especialmente si el n´ umero de grados de libertad es tambi´en elevado. Dentro de este grupo se encuentran los m´etodos denominados Mapa de Carreteras, Descomposici´on de Celdas y los de Campos de Potencial, seg´ un la clasificaci´on de Latombe [62]: Mapas de carretera.- Este m´etodo captura el espacio libre mediante una

8

´ CAP. 1. INTRODUCCION denominada red de carreteras a la cual se conectan los puntos inicial y final, obteni´endose la trayectoria como una ruta dentro de dicho mapa. Descomposici´on en celdas.- Se basa en la descomposici´on del espacio libre en zonas o celdas, de tal manera que una trayectoria entre un punto inicial y otro final se encontrar´a a trav´es de la sucesi´on de celdas contiguas que los une, siguiendo un m´etodo de b´ usqueda. Funciones potenciales y funciones potenciales arm´onicas.- Los m´etodos basados en funciones potenciales consideran el robot, representado como un punto en el espacio de configuraciones, como una part´ıcula de carga sometida a la fuerza de un campo de potencial creado artificialmente: los obst´aculos generar´an fuerzas repulsivas y el punto destino generar´a una fuerza atractiva; de esta forma la part´ıcula encontrar´a la trayectoria, libre de colisiones, siguiendo el gradiente del campo de potencial.

El m´etodo de las funciones potenciales, presentado por Khatib [56], tiene el inconveniente de que puede presentar m´ınimos locales, en donde la part´ıcula de carga que representa al robot quedar´ıa atrapada. Para intentar escapar de estos m´ınimos, el planificador, debe generar caminos aleatorios; tal como se expone en [62], y donde tambi´en se muestran los dos m´etodos num´ericos utilizados para generar las funciones de navegaci´on: Funci´on Num´erica Simple de Navegaci´on (NF1) y Funci´on Num´erica Mejorada de Navegaci´on (NF2). La funci´on NF1 se computa lanzando un frente de onda equipotencial desde el punto destino, tiene el inconveniente de generar trayectorias que pasan rozando a los obst´aculos, mientras que la funci´on NF2 se genera lanzando el frente de onda desde el esqueleto del C-espacio, computado previamente. Connolly, Burns and Weiss R. [25] plantean la utilizaci´on de la funci´on potencial como soluci´on de la ecuaci´on de Laplace. Las funciones que son soluci´on de la ecuaci´on de Laplace se llaman funciones arm´onicas, y el potencial basado en dichas funciones tiene la importante propiedad de no presentar m´ınimos locales. El c´alculo num´erico de la funci´on potencial arm´onica lo realiza Masoud en [73] utilizando un grafo, cuyos nodos representan la partici´on conexa del espacio de configuraciones, como si fuese una red el´ectrica y aplicando la segunda ley de Kirchoff (suma de corrientes igual a cero en los nodos). As´ı, el potencial en cualquier nodo toma el valor medio de todos sus vecinos, lo que implica la ausencia de m´ınimos locales para el potencial en los nodos. De forma similar lo hace Althofer, Fraser y Bugmann en [3]. Mientras que Garrido en [35] utiliza la teor´ıa de elementos finitos sobre una malla de c´elulas triangulares. Sobre los obst´aculos se aplican las condiciones de contorno (Neuman o Dirichlet) necesarias para la existencia de una soluci´on. Por ejemplo, en [74] a las regiones obst´aculo le adjudica valor cero (condiciones de contorno de Dirichlet). Tambi´en se aplican las funciones arm´onicas a entornos din´amicos, como en [75, 76] donde se presenta un control din´amico de un planificador de trayec-

´ CAP. 1. INTRODUCCION

9

torias basado en funciones arm´onicas, aplicado a un entorno din´amico cuyos obst´aculos son detectados mediante sensores. Se obtiene una se˜ nal de control del navegador sumando al gradiente del potencial un t´ermino proporcional a la velocidad de desplazamiento. Por su parte Hussein y Elnagax en [45], utilizan las ecuaciones de Maxwell para obtener una funci´on potencial que se adapte a entornos din´amicos y que no presente zonas donde el m´odulo de potencial presenta valores excesivamente peque˜ nos.

1.2.2.2.

M´ etodos basados en muestreo

Con la finalidad de evitar el c´alculo de los C-obst´aculos, y siguiendo procedimientos deterministas o aleatorios, estos m´etodos comprueban un n´ umero determinado de puntos del espacio de configuraciones; clasific´andolos como libres o de obst´aculo y obteni´endose una captura parcial del espacio de configuraciones. Con esta informaci´on, siguiendo distintos m´etodos, se busca una trayectoria; si existe una soluci´on la encuentran con una determinada probabilidad, que se acercar´a a uno seg´ un el n´ umero de exploraciones se acerca a infinito. Tienen la ventaja de ser aplicables para un n´ umero de dimensiones (grados de libertad) elevados. Dentro de este apartado se encuentran los m´etodos siguientes: Mapa de carreteras probabil´ıstico o Probabilistic Roadmap Method (PRM) Latombe y Kavraki [52].- Crean un mapa de caminos con puntos de cruce o nodos, los cuales se determinan muestreando aleatoriamente el espacio de configuraciones, utilizando para ello funciones detectoras de colisi´on. Se necesita un planificador local para enlazar los nodos y los puntos inicial y final deben de enlazarse previamente con el mapa. El m´etodo PRM es un planificador de movimientos muy utilizado, que tiene buen comportamiento en los robots con muchos grados de libertad. Van den Berg y Overmars en [8] plantean, sin embargo, que el m´etodo PRM no funciona tan bien en situaciones en las que el robot tiene que pasar a trav´es de pasillos largos y estrechos, en el espacio de configuraciones. Esto se debe principalmente a la uniformidad de las muestras utilizadas en el planificador, de manera que en las grandes regiones abiertas caen muchas muestras y muy pocas en los pasillos estrechos. Existen diferentes m´etodos para sesgar el muestreo. Por ejemplo Kazemi y Mehrandezh en [54], proponen mejorar la eficiencia del m´etodo PRM mediante la utilizaci´on de una funci´on potencial arm´onica; la cual se utiliza para condicionar la zona de muestreo, realizando m´as muestras en aquellas zonas se˜ naladas por el camino m´as corto obtenido por la funci´on arm´onica. De forma similar, Aamo, Kragic y Christensen en [1], proponen mejorar la eficiencia del m´etodo PRM, utilizando una funci´on potencial artificial para polarizar o condicionar el muestreo aleatorio y conducirlo hacia regiones de pasillos. Alternativamente, se ha planteado el uso de muestreo determinista en

10

´ CAP. 1. INTRODUCCION vez de muestreo aleatorio. Por ejemplo, Branicky et al, en [14], proponen el m´etodo cuasi aleatorio, Quasi-random Roadmap (Q-PRM), que mejora la dispersi´on de las muestras obtenidas en el PRM, mediante la secuencia determinista Hamniersley, utilizada en algunos m´etodos num´ericos. De esta manera consiguen solventar el problema persistente de los pasillos largos y estrechos, en espacios de un n´ umero de grados de libertad elevado. ´ Arboles de exploraci´on aleatoria r´apida o Rapidly Exploring Random Trees (RRT), LaValle [63].- Realiza una captura del espacio libre mediante una topolog´ıa en forma de ´arbol que evoluciona, con una resoluci´on creciente, utilizando t´ecnicas de muestreo aleatorio o determinista y funciones detectoras de colisi´on. Realiza una exploraci´on de puntos que, si son libres, va incorporando a la estructura mediante segmentos, utilizando un algoritmo que comprueba la proximidad al punto o segmento m´as cercano. En [64] se presenta el m´etodo b´asico, mientras que en [58] se propone la mejora de avanzar todo lo posible, desde la hoja m´as cercana, en la direcci´on de la muestra sacada aleatoriamente. En [100] se limita el espacio donde se muestrea, para garantizar que el crecimiento del ´arbol es posible. Descomposici´on aleatoria de celdas o Probabilistic Cell Decomposition (PCD), Lingelbach [67].- Realizan una descomposici´on del espacio en celdas, las cuales son muetreadas aleatoriamente y descompuestas en 2n -tree, si se encuentran puntos libres y de obst´aculo. Sobre estas celdas se realizan b´ usquedas en ´arbol hasta que se encuentra una trayectoria, o bien ha transcurrido un tiempo determinado. Tambi´en Lindemann y LaValle en [66], proponen la utilizaci´on de una secuencia determinista de muestreo que presenta mejores propiedades de cobertura uniforme que la aleatoria. Con dicho muestreo se obtiene una estructura de rejillas en multirresolici´on, la cual permite una adaptaci´on ´optima de la descomposici´on en celdillas. Burlet, Aycard y Fraichard en [18], proponen la combinaci´on de la descomposici´on de celdas en quad-tree con la generaci´on de cadenas de Markov para la obtenci´on de trayectorias. Cada celda que forma parte de la cadena se muestrea y, si contiene obst´acuo, se descompone y se repite el procedimiento. Combinaci´on de m´etodos aleatorios y funciones potenciales.- Una de las primeras aplicaciones de los m´etodos aleatorios, en la planificaci´on de movimientos en rob´otica, fue la generaci´on de caminos aleatorios para intentar escapar de los m´ınimos locales presentes en las funciones potenciales, utilizadas en los planificadores basados en las mismas [62]. Otro enfoque es utilizar una funci´on potencial para condicionar el muestreo aleatorio en los m´etodos PRM, favoreciendo la exploraci´on de las regiones con pasillos [1], tal como se ha comentado anteriormente.

´ CAP. 1. INTRODUCCION 1.2.2.3.

11

M´ etodos para entornos din´ amicos

Estos m´etodos tienen en cuenta las variaciones temporales del entorno de trabajo, bien por el movimiento de los obst´aculos o por el cambio de escenario. Teniendo en cuenta la velocidad de cambio del entorno en relaci´on con el tiempo de ejecuci´on de las trayectorias, existen dos tipos de m´etodos para entornos din´amicos: los que consideran la posibilidad de colisi´on y los que no tienen en cuenta la posibilidad de colisi´on: Los planificadores que tienen en cuenta la posibilidad de colisi´on consideran que los desplazamiento de los obst´aculos pueden intersectar con la trayectoria del robot. No tienen conocimiento de la din´amica del entorno, disponen de un sistema sensorial y se basan en un sistema de control de movimientos reactivos [26]. Los planificadores para entornos din´amicos que no tienen en cuenta la posibilidad de colisi´on se basan en la predicci´on, con un cierto grado de certidumbre, de las posiciones futuras de los obst´aculos en movimiento, para lo cual precisan de una detecci´on de las posiciones y velocidades de los mismos, as´ı como de la memorizaci´on de datos hist´oricos. Con este enfoque, Gao y Sun [34] proponen un planificador para robots m´oviles desplaz´andose en entornos din´amicos, realizando previsiones de los movimientos de los obst´aculos y definiendo zonas peligrosas de colisi´on mediante la utilizaci´on de un filtro Kalman. Por su parte, Metoui et al [78] proponen un planificador robusto para robots m´oviles en un entorno de obst´aculos din´amicos. Este planificador considera una funci´on potencial con una fuerza de atracci´on al punto destino, fuerzas de repulsi´on de los obst´aculos y una fuerza atractiva del robot m´ovil respecto de los obst´aculos que tiene en cuenta tanto la posici´on como velocidad relativa.

1.3.

Objetivos

La utilizaci´on de las funciones arm´onicas, para generar funciones potenciales, evita la aparici´on de m´ınimos locales y estas se adaptan bien a todo tipo de entornos: angostos, din´amicos e inciertos. Sin embargo, este m´etodo usa habitualmente una discretizaci´on regular del espacio de configuraciones en forma de rejilla; sobre ella se realiza el c´alculo num´erico de la funci´on y su paso debe ser el de la precisi´on especificada. De esta forma, si la precisi´on y el n´ umero de grados de libertad son elevados, el m´etodo es inviable por el tiempo de computaci´on requerido. El objetivo de esta tesis es superar esta limitaci´on, dise˜ nando m´etodos y algoritmos que consigan aplicar las funciones arm´onicas en la realizaci´on de un Planificador de movimientos de manera m´as eficiente, permitiendo abordar problemas con un n´ umero de grados de libertad m´as elevado. A pesar del problema de eficiencia de computaci´on de las funciones arm´onicas, su c´alculo num´erico se reduce a una suma ponderada con un cierto n´ umero de iteraciones, lo que supone un coste computacional m´ınimo. Adem´as, el caso de entornos complicados plantea una dificultad que es com´ un a todos los m´etodos, la obtenci´on de trayectorias atrav´es de pasillos largos y estrechos.

12

´ CAP. 1. INTRODUCCION

Figura 1.2: Organigrama para la consecuci´on de los objetivos. Por ello, con este trabajo, se pretende mejorar la metodolog´ıa para hacer factible el uso de las funciones arm´onicas en el desarrollo de un planificador de movimientos. As´ı, de acuerdo con la figura 1.2, se proponen dos objetivos fundamentales con una serie de subobjetivos: Objetivo 1.- Estudio te´orico de las funciones arm´onicas como herramienta para la planificaci´on de movimientos: Estudio de las propiedades de las funciones arm´onicas. Propuesta de moldeado para evitar inconvenientes. Evaluaci´on de la aplicabilidad en entornos din´amicos. Objetivo 2.- Estudio de m´etodos para la implementaci´on eficiente de planificadores basados en funciones arm´onicas. Estudio de m´etodos num´ericos y propuesta para el c´alculo eficiente de funciones arm´onicas de cara a la planificaci´on de movimientos. Estudio de m´etodos basados en muestreo para la exploraci´on del C-espacio y su posible combinaci´on con m´etodos basados en funciones arm´onicas. El c´alculo de trayectorias requiere m´etodos que reduzcan dr´asticamente el n´ umero de puntos de c´alculo sin que por ello se pierda precisi´on. La utilizaci´on de las funciones arm´onicas en rejillas no uniformes, que discretizan el espacio de configuraciones de forma recursiva y jerarquizada, es ´optima para este fin. Este m´etodo de discretizaci´on se basa en la partici´on 2n -tree que de forma iterativa se adapta y define los contornos de los obst´aculos, dejando los espacios libres y de obst´aculos con celdas gruesas mientras los contornos se definen con celdas de paso fino, requerido por la precisi´on. Sin embargo, esta

´ CAP. 1. INTRODUCCION

13

discretizaci´on no se utilizar´a para realizar una b´ usqueda si no para realizar sobre ´el el c´alculo de las funciones arm´onicas. En previsi´on de imprecisiones y posibles variaciones en la posici´on del obst´aculo, es necesario desarrollar una metodolog´ıa que permita una adaptaci´on en l´ınea, la cual es posible gracias a la rapidez de c´alculo de actualizaci´on del valor de la funci´on arm´onica, cuando cambian ligeramente las condiciones de contorno. Por otro lado los m´etodos aleatorios, aplicados a la planificaci´on de movimientos, permiten encontrar una ruta sin necesidad de realizar una exploraci´on exhaustiva del espacio de configuraciones. En este ´ambito, los m´etodos aleatorios con sesgo de exploraci´on han mostrado ser m´as eficientes que los puramente aleatorios, ya que en ellos existe una tendencia de base que permite un mejor comportamiento en canales estrechos y largos. As´ı, los m´etodos aleatorios condicionados o guiados por el campo de potencial, pueden realizar una exploraci´on preferente de las zonas m´as prometedoras, acelerando el proceso. La utilizaci´on coordinada de los distintos m´etodos comentados, se presenta como un procedimiento id´oneo para desarrollar un planificador de movimientos u ´til y eficiente.

14

´ CAP. 1. INTRODUCCION

Cap´ıtulo 2 Planificaci´ on con funciones potenciales arm´ onicas, subarm´ onicas y superarm´ onicas La planificaci´on de movimientos de un robot se realiza usualmente en el espacio de configuraciones (C-espacio), donde una configuraci´on del robot se representa por un punto q, pudiendo corresponder a una posici´on libre o de colisi´on con un determinado obst´aculo; de esta forma los obst´aculos del espacio f´ısico son representados como obst´aculos en el espacio de configuraciones (C-obst´aculos). Tambi´en cabe la posibilidad de que el robot colisione consigo mismo, en este caso la zona de colisi´on sobre el espacio de configuraciones tiene la misma consideraci´on que los dem´as obst´aculos. As´ı, el movimiento del robot entre dos configuraciones, libre de colisiones con los obst´aculos, se traduce en una trayectoria curvil´ınea del punto q sin intersecciones con los C-obst´aculos. La planificaci´on mediante funciones potenciales consiste en definir una funci´on sobre el espacio libre de configuraciones, cuyo gradiente act´ ue como un campo potencial que arrastre cualquier carga puntual situada en su dominio hasta un m´ınimo global donde se sit´ ua el punto destino. De esta forma, la trayectoria se obtiene mediante el seguimiento del gradiente, sobre la hipersuperficie de la funci´on potencial, entre el punto inicial y final. La funci´on potencial creada como suma de una repulsiva respecto a los obst´aculos y otra de atracci´on hacia el objetivo, tiene el problema de presentar m´ınimos locales [56], mientras que la funci´on potencial creada como soluci´on de la ecuaci´on de Laplace [25] no presenta m´ınimos locales y es la base del presente trabajo. En este cap´ıtulo se analizan las propiedades de las funciones arm´onicas, subarm´onicas y superarm´onicas, desde el punto de vista de su aplicaci´on en la planificaci´on de movimientos en rob´otica. Se destacan tanto aquellas propiedades que son favorables, como las que son desfavorables y han limitado su aplicaci´on a robots con un n´ umero de grados de libertad reducido. Con el mismo enfoque se analizan tambi´en las funciones subarm´onicas y superarm´onicas, averiguando sus posibilidades en nuevas aplicaciones para la planificaci´on de movimientos. 15

´ CAP. 2. FUNCIONES ARMONICAS

16

Figura 2.1: Planificaci´on en el espacio de configuraciones, considerando los obst´aculos como parte del contorno.

2.1.

Introducci´ on

Una funci´on potencial es una funci´on escalar u(r) definida en un dominio Ω ⊂ Rn . Para la planificaci´on de trayectorias Ω es coincidente con el espacio libre de configuraciones y el gradiente de dicha funci´on constituye un campo vectorial ∇u = (∂u/∂q1 , . . . , ∂u/∂qn ) (2.1) que, para que sea u ´til en la planificaci´on de trayectorias, no debe tener puntos sumidero de flujo del gradiente excepto en el punto destino PD situado en rd . Es decir, la funci´on potencial u debe cumplir que div (grad u) ≤ 0,

∀r ∈ (Ω − PD )

(2.2)

Lo que permite dos posibilidades en el tratamiento o modelado de los obst´aculos: I. Considerar los obst´ aculos junto con las condiciones de contorno. El dominio Ω no contiene los obst´aculos y una funci´on delta de Dirac se sit´ ua en el punto destino (figura 2.1). El contorno del espacio libre Ω se obtiene como la uni´on del contorno externo y los contornos de los obst´aculos ∂Ω = ∂ΩE ∪ ∂O1 ∪ ∂O2 ∪ ∂O3 . . .

(2.3)

La funci´on potencial se obtiene como soluci´on de la ecuaci´on de Poisson ∇2 u = δ(r − rd ),

(2.4)

de forma que las l´ıneas de corriente definidas por el gradiente de la funci´on u, siempre terminar´an en el punto situado en rd . Es decir, la funci´on potencial

´ CAP. 2. FUNCIONES ARMONICAS

17

Figura 2.2: Planificaci´on en el espacio de configuraciones, considerando los obst´aculos como una funci´on. solo presentar´a un m´ınimo, en todo el espacio de configuraciones, situado en el punto de aplicaci´on de la funci´on delta. Por otro lado, para cualquier subdominio que no contenga el punto rd se cumplir´a la ecuaci´on de Laplace ∇2 u = 0,

(2.5)

cuya soluci´on es una funci´on arm´onica que tiene la importante propiedad de no presentar m´ınimos locales. Es por esto, que la construcci´on de una funci´on potencial arm´onica servir´a para obtener un algoritmo de navegaci´on que, mediante la b´ usqueda del m´ınimo global, reporte la trayectoria deseada.

II. Considerar los obst´ aculos como una funci´ on ρ. En este caso el dominio Ω contiene los obst´aculos, igualmente se sit´ ua una funci´on delta en el punto destino (figura 2.2) y el subdominio obst´aculo Ωo se compone de la uni´on de los subdominios de los obst´aculos Ωo = Ωo1 ∪ Ωo2 ∪ Ωo3 . . .

(2.6)

As´ı, el modelo planteado es ∇2 u = −ρ(r) + δ(r − rd ), la funci´on ρ representa los obst´aculos de forma que   > 0 ∀r ⊂ ΩO ρ(r) =  = 0 ∀r * ΩO

(2.7)

´ CAP. 2. FUNCIONES ARMONICAS

18

La soluci´on de (2.7) ser´a una funci´on superarm´onica que contendr´a un n´ umero determinado de m´aximos y ning´ un m´ınimo. Por consiguente, en todo el dominio solo existir´a un m´ınimo situado en rd . En cualquiera de los dos casos un algoritmo auxiliar de navegaci´on obtendr´a la trayectoria siguiendo el gradiente hasta alcanzar el objetivo. Esta trayectoria consistir´a en una curvilinea (unidimensional), coincidente con una l´ınea de corriente, la cual se puede representar de forma param´etrica lc = r(s),

s ∈ [si , sf ]

(2.8)

donde si y sf son los valores del par´ametro en el punto inicial y final. El vector tangente, en cada punto, ser´a igual al gradiente de la funci´on t(s) =

dr = ∇u ds

(2.9)

y cuya curvatura, tambi´en en cada punto, ser´a k(s) =

2.2.

dt d(∇u) = ds ds

(2.10)

Propiedades afines a la generaci´ on de funciones potenciales para la planificaci´ on de movimientos en rob´ otica

Los entornos de trabajo en rob´otica son complicados, de forma que la soluci´on de la ecuaci´on de Laplace o de Poisson se obtendr´a por m´etodos num´ericos; siendo posible la obtenci´on de una expresi´on anal´ıtica de u solamente para modelos sencillos. En cualquier caso lo que importa, en este punto del estudio, es conocer las propiedades de las soluciones m´as que las propias soluciones, ya que as´ı se puede conocer la viabilidad y condicionantes de la aplicaci´on de las funciones arm´onicas a la planificaci´on de movimientos. Los aspectos y propiedades principales a analizar de las funciones arm´onicas, aplicables al desarrollo de un planificador de trayectorias, son: Ausencia de m´ınimos locales. Concavidad, convexidad de las funciones arm´onicas y superarm´onicas. L´ıneas de corriente y superficies equipotenciales. Regularidad de la soluci´on; curvatura de las trayectorias. Condiciones de contorno, condici´on de solubilidad y balance global del gradiente. Estabilidad de la soluci´on. Superposici´on lineal.

´ CAP. 2. FUNCIONES ARMONICAS

19

Figura 2.3: Principio de localidad de una funci´on de dos dimensiones.

2.2.1.

Ausencia de m´ınimos locales

Esta propiedad es esencial para la utilizaci´on de las funciones arm´onicas y superarm´onicas en la generaci´on de funciones potenciales, utilizadas en la planificaci´on de trayectorias. Se cumple si y solo si para cualquier entorno pr´oximo a todo punto del dominio, el valor m´ınimo de la funci´on se encuentra sobre su frontera. Los entornos considerados ser´an hiperesferas centradas en el punto, cuyo radio podr´a ser tan peque˜ no como se quiera; sobre estos entornos se analizan las propiedades del valor medio, del m´aximo y del m´ınimo. 2.2.1.1.

An´ alisis en dos dimensiones

Como paso previo a la generalizaci´on para n dimensiones, se realiza el estudio en dos dimensiones, de esta forma se puede ilustrar los resultados mostrando una visi´on m´as intuitiva de los mismos. Proposici´ on 1. Principio de localidad de una funci´ on: Dada una 2 2 funci´on f (x, y), definida en un dominio Ω ⊂ R y de clase C . Dado un disco de radio R centrado en (x0 , y0 ): D = {(x, y) ∈ Ω | k(x − x0 , y − y0 )k < R}. Entonces se cumple que: I R2 1 div (grad f )0 (2.11) f (x, y)dl − f (x0 , y0 ) = 2πR L 4 donde la integral de l´ınea se realiza sobre la circunferencia contorno del disco, tal como se muestra en la figura 2.3. Esta sugerente expresi´on dice que el valor de una funci´on (de clase C 2 ) en un punto es igual a su valor medio sobre una circunferencia centrada en dicho punto de radio R, menos un t´ermino proporcional a la divergencia del gradiente (laplaciana) en el mismo punto.

´ CAP. 2. FUNCIONES ARMONICAS

20

Figura 2.4: Entorno de un punto (x0 , y0 ) en una funci´on arm´onica. La demostraci´on de esta proposici´on se realiza en el ap´endice A.1. Funciones arm´ onicas: En el caso de funciones arm´onicas se cumple la ecuaci´on de Laplace: ∂ 2f ∂ 2f + =0 ∂x2 ∂y 2

⇒

div (grad f ) = 0 ∀(x, y) ∈ Ω

(2.12)

Siendo fM ED(L) , fM IN (L) y fM AX(L) los valores medio, m´ınimo y m´aximo de f sobre la circunferencia de radio R, respectivamente, la ecuaci´on (2.11) da I I fM ED(L) 1 f (x0 , y0 ) = f (x, y)dl = dl = fM ED(L) (2.13) 2πR L 2πR L

Cumpli´endose la desigualdad I I I 1 1 1 f (x, y)dl < fM IN (L) dl < fM AX(L) dl 2πR 2πR L 2πR L L fM IN (L) < f0 < fM AX(L)

(2.14) (2.15)

Es decir, el valor de la funci´on en el centro es el valor medio de los valores sobre la circunferencia y, por tanto, esta comprendido entre el m´ınimo y el m´aximo que se encuentran sobre la circunferencia (figura 2.4). De esta forma se puede formar una sucesi´on concatenada de esferas, tal que cada una contenga en su interior el valor m´ınimo de la esfera precedente, de tal forma que esta sucesi´on de esferas (figura 2.5) termina cuando la u ´ltima tiene el m´ınimo sobre el contorno del dominio, lo que supone la ausencia de m´aximos y m´ınimos locales en todo el dominio. Por ello, las trayectorias generadas mediante el seguimiento del gradiente de una funci´on potencial arm´onica, definida en el dominio del espacio de configuraciones con un contorno ∂Ω y una funci´on δ situada en el punto objetivo, no caer´an en m´ınimos locales y tendr´an siempre soluci´on. Cabe la posibilidad de existencia de puntos de ensilladura, pero estos no suponen un problema porque el seguimiento del gradiente siempre encontrar´ıa una salida, siguiendo una de las vertientes. Funciones subarm´ onicas: Las funciones subarm´onicas se definen como aquellas que son soluci´on de la ecuaci´on de Poisson ∂ 2f ∂ 2f + = ρ(x, y), ∂x2 ∂y 2

(2.16)

´ CAP. 2. FUNCIONES ARMONICAS

21

Figura 2.5: M´ınimo en el contorno para una funci´on arm´onica. tal que ρ(x, y) ≥ 0 ∀(x, y) ∈ Ω

(2.17)

div (grad f ) = ρ(x, y) ≥ 0 ∀(x, y) ∈ Ω

(2.18)

Es decir Sustituyendo esta condici´on en la expresi´on del Principio de localidad de una funci´on (ecuaci´on 2.11), se obtiene I 1 1 f (x0 , y0 ) = ρ(x0 , y0 )R2 (2.19) f (x, y)dl − 2πR L 2 · 2! de donde se deduce la desigualdad 1 fM IN (L) 2πR

I

L

1 ρ(x0 , y0 )R2 < 2 · 2! I 1 1 ρ(x0 , y0 )R2 < f (x, y)dl − 2πR L 2 · 2! I 1 1 fM AX(L) dl − ρ(x0 , y0 )R2 (2.20) 2πR 2 · 2! L

dl −

fM IN (L) − λ < f0 < fM AX(L) − λ,

λ=

1 ρ(x0 , y0 )R2 ≥ 0 2 · 2!

(2.21)

Es decir, el valor m´aximo de la funci´on se encuentra en la frontera y el m´ınimo en el interior o en la frontera (figura 2.6). Funciones superarm´ onicas: Las funciones subarm´onicas son soluci´on de ∂ 2f ∂ 2f + = −ρ(x, y) | ρ(x, y) ≥ 0 ∀(x, y) ∈ Ω (2.22) ∂x2 ∂y 2 Y sustituyendo div (grad f ) = −ρ(x, y)

(2.23)

´ CAP. 2. FUNCIONES ARMONICAS

22

Figura 2.6: Funciones subarm´onicas y superarm´onicas en la ecuaci´on 2.12, se tiene 1 f (x0 , y0 ) = 2πR

I

f (x, y)dl + l

1 ρ(x0 , y0 )R2 2 · 2!

(2.24)

de donde se deduce I 1 1 fM IN (L) dl + ρ(x0 , y0 )R2 < 2πR 2 · 2! L I 1 1 ρ(x0 , y0 )R2 < f (x, y)dl + 2πR L 2 · 2! I 1 1 fM AX(L) dl + ρ(x0 , y0 )R2 (2.25) 2πR 2 · 2! L fM IN (L) + λ < f0 < fM AX(L) + λ,

λ=

1 ρ(x0 , y0 )R2 ≥ 0 2 · 2!

(2.26)

Es decir el valor m´ınimo de la funci´on se encuentra en la frontera y m´aximo en el interior o en la frontera (figura 2.6) Esta importante propiedad permite la utilizaci´on de las funciones superarm´onicas, igual que las arm´onicas, para la generaci´on de trayectorias libres de colisiones y sin m´ınimos locales.

2.2.1.2.

An´ alisis en n dimensiones

Considerar una funci´on potencial de n variables: u = f (q1 , q2 , . . . , qn ) ,

(2.27)

anal´ıtica en el dominio Ω ⊂ Rn , con un contorno ∂Ω sobre el que se impondr´an unas condiciones a la funci´on. Las coordenadas se pueden expresar vectorialmente, r = (q1 , q2 , . . . , qn ) y la distancia euclidiana entre dos puntos como d = kr1 − r2 k

´ CAP. 2. FUNCIONES ARMONICAS

23

Generalizaci´ on del Principio de localidad de una funci´ on: Dada una funci´on f (r), definida en un dominio Ω ⊂ Rn y de clase C ∞ . Dado un disco de radio R centrado en (r0 ): D = {(r) ∈ Ω | kr − r0 k < R}. Entonces se cumple que: 1 f (r0 ) = S

I

S

f (r0 + Rn)dS −

∞ X n=1


R2n ∇2n f r0 2 · (2n)!

(2.28)

Esta expresi´on es la generalizaci´on de la 2.11 y la demostraci´on de su validez se muestra en el ap´endice A.1. Funciones arm´ onicas: Las funciones arm´onicas cumplen ∇2 f = 0, en todo punto perteneciente al dominio, aplicando esta condici´on al principio de localidad (ecuaci´on 2.28), se obtiene la propiedad del valor medio para una funci´on multivariada: I 1 f (r0 + Rn)dS (2.29) f (r0 ) = S S

Por lo tanto, el valor de una funci´on arm´onica en cualquier punto de un volumen determinado no ser´a ni m´aximo ni m´ınimo, encontr´andose estos sobre la superficie cerrada que limita dicho volumen: H H H fM IN (S) S dS f dS f (r + Rn)dS 0 M AX(S) S < S < (2.30) S S S fM IN (S) < f (r0 ) < fM AX(S)

(2.31)

Funciones subarm´ onicas: Igual que para el caso de dos dimensiones, las funciones subarm´onicas para n dimensiones son soluci´on de ∇2 f = ρ (q1 , q2 , . . . , qn )

| ρ (q1 , q2 , . . . , qn ) ≥ 0 ∀ (q1 , q2 , . . . , qn ) ∈ Ω (2.32) Aplicado al principio de localidad (ecuaci´on 2.28), se tiene: H f (r0 + Rn)dS R2 −ρ (2.33) f (r0 ) = S S 4 Con respecto a los valores m´aximo y m´ınimo en la hiperesfera, se cumple la desigualdad: H H H f (r0 + Rn)dS fM IN (S) S dS fM AX(S) S dS R2 R2 R2 S −ρ < −ρ < −ρ S 4 S 4 S 4 (2.34) R2 R2 < f (r0 ) < fM AX(S) − ρ (2.35) fM IN (S) − ρ 4 4 Es decir, m´aximo en la frontera y m´ınimo en el interior o en la frontera (figura 2.6)

´ CAP. 2. FUNCIONES ARMONICAS

24

Figura 2.7: Funci´on convexa.

r0 = r1 (1 − ϑ) + r2 ϑ,

0≤ϑ≤1

Funciones superarm´ onicas: Similarmente, las funciones superarm´onicas cumplen ∇2 f = −ρ (q1 , q2 , . . . , qn )

| ρ (q1 , q2 , . . . , qn ) ≥ 0 ∀ (q1 , q2 , . . . , qn ) ∈ Ω (2.36) H 2 f (r + Rn)dS R 0 +ρ (2.37) f (r0 ) = S S 4

Y la desigualdad que involucra a los valores m´aximo y m´ınimo H H H fM IN (S) S dS fM AX(S) S dS f (r0 + Rn)dS R2 R2 R2 S +ρ < +ρ < +ρ S 4 S 4 S 4 (2.38) V V (2.39) fM IN (S) + ρR < f (r0 ) < fM AX(S) + ρR S S De donde se deduce que el m´ınimo se sit´ ua en la frontera y el m´aximo en el interior o en la frontera (figura 2.6).

2.2.2.

Concavidad, convexidad de las funciones subarm´ onicas y superarm´ onicas

La propiedad del valor medio y del m´ınimo, que cumplen las funciones arm´onicas y superarm´onicas, posibilita su utilizaci´on en la generaci´on de funciones potenciales definidas en el dominio del espacio de configuraciones. Para saber m´as respecto al comportamiento del campo vectorial determinado por el gradiente de dicha funci´on vectorial es necesario estudiar la concavidad y convexidad de dicha funci´on. As´ı, el cumplimiento de la desigualdad de Jensen [13] (ecuaci´on 2.41 y figura 2.7) determina la no concavidad de las funciones arm´onicas y superarm´onicas, lo cual significa que la variaci´on

´ CAP. 2. FUNCIONES ARMONICAS

25

del m´odulo del gradiente no cambiar´a de decreciente a creciente en ninguna direcci´on. Por consiguiente tampoco existir´an m´ınimos locales en todo el dominio de la funci´on. P 2 Por lo tanto, para ni=1 ∂∂xf2 ≥ 0, se cumple la desigualdad de Jensen para i las funciones c´oncavas: f (r1 (1 − ϑ) + r2 ϑ) ≤ f (r1 )(1 − ϑ) + f (r2 )ϑ,

0≤ϑ≤1

(2.40)

P 2 Y para ni=1 ∂∂xf2 ≤ 0, se cumple la desigualdad de Jensen para las funciones i convexas (figura 2.7): f (r1 (1 − ϑ) + r2 ϑ) ≥ f (r1 )(1 − ϑ) + f (r2 )ϑ,

0≤ϑ≤1

(2.41)

donde r1 y r2 son dos puntos cualesquiera del espacio. La demostraci´on de esta propiedad se muestra en el ap´endice A.3.1.

2.2.3.

L´ıneas de corriente y superficies equipotenciales

En todos los puntos de una superficie equipotencial Seq , el valor de la funci´on es el mismo: f (r) = Const ∀r ∈ Seq . Por lo que la derivada direccional, en todo punto de Seq , con respecto a cualquier vector v contenido en dicha superficie ser´a cero: ∂f = ∇f · v = 0 ∂v

→

∇f ⊥v

(2.42)

por lo que las l´ıneas de corriente y las superficies equipotenciales son tambi´en ortogonales. Adem´as la ecuaci´on de Laplace dice que ∇2 u = div (grad u) = 0,

(2.43)

lo que determina que las l´ıneas de corriente no se cortar´an, ya que la divergencia del gradiente es nula en todo punto del dominio. Para un entorno del espacio de configuraciones modelado por la ecuaci´on de Poisson, con condiciones de contorno definidas sobre los obst´aculos y una funci´on delta independiente situada en el punto objetivo, la funci´on soluci´on dispondr´a de infinitas l´ıneas de corriente que confluir´an en el punto objetivo. De esta forma la trayectoria, libre de colisiones, se superpondr´a sobre una de estas l´ıneas que, partiendo del punto configuraci´on inicial, llegar´a a dicho punto objetivo. La posibilidad de corregir la trayectoria pasando de una l´ınea de corriente a otra pr´oxima, a trav´es de una superficie equipotencial, ofrece un m´etodo para dise˜ nar un control reactivo para abordar entornos din´amicos o inciertos.

´ CAP. 2. FUNCIONES ARMONICAS

26

Figura 2.8: Estudio de la curvatura.

2.2.4.

Superposici´ on lineal

Por la propiedad de linealidad del laplaciano se deduce que la combinaci´on lineal de varias funciones arm´onicas es, a su vez, una funci´on arm´onica: Dado un dominio Ω ⊂ Rn y dada una serie de funciones ui : Ω → R, donde i = 1, 2, 3, . . . , k. Si cada una de las funciones ui es arm´onica en el dominio, entonces la combinaci´on lineal de ellas tambi´en lo es ! k k X X
2 λ i ∇ 2 ui = 0 (2.44) λ i ui = ∇ i=1

i=1

Esta propiedad permite aplicar el principio de superposici´on a la soluci´on de la ecuaci´on de Laplace o de Poisson; tal como se hace con las funciones de Green, lo que facilita el estudio de la influencia de las condiciones de contorno, as´ı como del modelado de los obst´aculos como funciones independientes en las ecuaciones sub y superarm´onicas.

2.2.5.

Curvatura de las funciones arm´ onicas

Una trayectoria sobre el espacio de configuraciones quedar´a definida por una curvil´ınea tal que el vector tangente en cada punto de la misma tendr´a la misma direcci´on que el gradiente de la funci´on en el mismo punto. La curvatura de esta linea, en el espacio de configuraciones, determinar´a las aceleraciones que experimentar´a el robot para la consecuci´on de la trayectoria. Interesa saber como son de “suaves”las trayectorias obtenidas a partir de las funciones arm´onicas. Resulta que para el caso de dichas funciones potenciales arm´onicas esta curvatura ser´a reducida, lo que proporcionar´a trayectorias con aceleraciones menores para los actuadores que impulsan las articulaciones, lo cual puede redundar en un menor consumo de energ´ıa y menos problemas de mantenimiento.

´ CAP. 2. FUNCIONES ARMONICAS

27

De acuerdo con la figura 2.8, considerando dos puntos pr´oximos de una curvilinea obtenida sobre una hipersuperficie de una funci´on potencial gen´erica, la curvatura de la curvil´ınea en la direcci´on v es    1 ∂f ∂f ∂ 2f K1v = l´ım (r2 ) − (r1 ) = (r1 ) (2.45) r2 →r1 |r2 − r1 | ∂v ∂v ∂v2 Entonces, teniendo en cuenta que la Laplaciana de una funci´on es la divergencia de su gradiente, y que por consiguiente ´esta se define en el entorno de un punto, se integra sobre una superficie esf´erica de radio muy peque˜ n o δR y se obtiene la curvatura media de Gauss:  I  1 1 K1M G = l´ım K1v dδS = ∇2 f, (2.46) δR −→0 δS δS 2 tal como se demuestra en el ap´endice A.3.2. Esta curvatura media ser´a m´ınima (cero) si se cumple la ecuaci´on de Laplace. Por lo tanto, utilizando las funciones arm´onicas para generar funciones potenciales auxiliares, obtendremos trayectorias “suaves”, lo que supone aceleraciones suaves para el robot en la ejecuci´on de movimientos.

2.2.6.

Estabilidad de la soluci´ on

La soluci´on de la ecuaci´on de Laplace o de Poisson experimenta peque˜ nas variaciones frente a peque˜ nas variaciones de las condiciones de contorno o peque˜ nas variaciones de la funci´on independiente, esta propiedad facilita la utilizaci´on de las funciones potenciales arm´onicas en el dise˜ no de planificadores de trayectorias en entornos cambiantes. En el ap´endice A.3.6 se analiza esta propiedad.

2.3.

Problem´ atica inherente al m´ etodo

Una vez obtenida la funci´on potencial, la trayectoria se obtiene por seguimiento del gradiente; por lo que la distribuci´on del mismo, y la existencia de zonas donde su m´odulo pueda alcanzar valores extremadamente bajos, puede presentar serios problemas para la computaci´on de trayectorias. Este efecto es caracter´ıstico de las funciones arm´onicas, est´a presente en la naturaleza (campos el´ectricos y gravitatorios) y supone la principal desventaja en su utilizaci´on para la planificaci´on de trayectorias en rob´otica. Con la finalidad de conservar una visi´on intuitiva del mismo se le denomina efecto embudo.

2.3.1.

Efecto embudo y funciones arm´ onicas

Con la finalidad de ilustrar el problema fundamental, presentado en la utilizaci´on de las funciones arm´onicas como funciones potenciales para la

´ CAP. 2. FUNCIONES ARMONICAS

28

Figura 2.9: Efecto embudo. generaci´on de movimientos en rob´otica, se considera el modelo simplificado siguiente: ∇2 u = δ (r − rd ) , (2.47) considerando todo el espacio y situando la funci´on delta en el origen de coordenadas (rd = 0). Aplicando el teorema de la divergencia al gradiente de la funci´on; para un hipervolumen esf´erico V centrado tambi´en en el origen de coordenadas, limitado por la hipersuperficie esf´erica correspondiente S, de radio R I Z Z ∇u · dS = (∇ · ∇u) dV = δ (r) dV = 1 (2.48) S

V

V

Por la simetr´ıa del problema ∇u · dS = |∇u|dS y |∇u| = constante, por lo que se tiene: I I ∇u · dS = |∇u| dS = |∇u|S, (2.49) S

S

sustituyendo en 2.48 se tiene que |∇u|S = 1 y ∇u =

1 ru Sn

(2.50)

Siendo ru el vector unitario en la direcci´on del radio de la hiperesfera, apuntando al exterior y aplicado a un punto de la superficie de la misma. La expresi´on de la superficie de una hiperesfera en n dimensiones se obtiene en el ap´endice A.2 y vale: n

π 2 nrn−1
Sn = Γ n2 + 1

(2.51)

As´ı, el m´odulo del gradiente decrece con el aumento del radio de la hiperesfera y este decrecimiento es m´as r´apido cuanto mayor es el n´ umero de variables de la funci´on potencial.

´ CAP. 2. FUNCIONES ARMONICAS

29

Figura 2.10: Balance del flujo del gradiente. De esta forma, para un n´ umero de grados de libertad elevado, el gradiente se hace tan peque˜ no que, a una cierta distancia r, la propia precisi´on num´erica del ordenador puede no ser suficiente para discriminarlo. Adem´as, el efecto embudo puede provocar que los m´etodos num´ericos de relajaci´on, utilizados en el c´alculo de las funciones arm´onicas, requieran un n´ umero muy elevado de iteraciones y un tiempo de computaci´on tambi´en muy elevado.

2.3.2.

Condiciones de contorno, condiciones de solubilidad y balance global del flujo del gradiente

Las condiciones de contorno aplicables al modelo, para la obtenci´on de la funci´on potencial, no pueden ser cualesquiera si no que deben cumplir con unas determinadas condiciones denominadas condiciones de solubilidad o de compatibilidad; las cuales surgen de aplicar el teorema de la divergencia al campo vectorial creado por el gradiente de la funci´on, en todo el dominio. Considerando el modelo gen´erico (ecuaci´on 2.7) ∇2 u = −ρ(r) + δ(r − rd ),

| ρ(r) > 0 ∀r ∈ Ω

(2.52)

aplicable a un dominio Ω con un contorno tambi´en gen´erico ∂Ω (figura 2.10). Considerando el campo vectorial g(r) = ∇u(r),

∀r ∈ Ω

aplicando a g el teorema de la divergencia Z I (∇ · g) dV = g · dS, V

(2.54)

S

considerando todo el dominio Ω y su contorno ∂Ω, se tiene Z I
2 ∇ u dΩ = ∇u · d(∂Ω), Ω

(2.53)

∂Ω

(2.55)

´ CAP. 2. FUNCIONES ARMONICAS

30

ahora, teniendo en cuenta que ∇u · n = Z

Ω

∂u ∂n

y 2.52

[−ρ(r) + δ(r − rd )] dΩ =

I

∂Ω

∂u d∂Ω, ∂n

(2.56)

de donde, se deduce I

∂Ω

∂u dS = 1 − ∂n

Z

ρ(r)dΩO

(2.57)

Ωo

que son las condiciones de solubilidad de la ecuaci´on 2.52 e imponen restricciones a las condiciones de contorno, que deben cumplirlas para que exista una soluci´on. Similarmente, aplicando 2.55 a un contorno formado por SC ∪ SO ∪ Sd , tal como se muestra en la figura 2.10, se tiene I I I ∇u · dSC − ∇u · dSO − ∇u · dSd = 0, (2.58) SC

SO

Sd

donde cada t´ermino representa los flujos de gradiente que atraviesa cada superficie φC = φO + φd , (2.59) expresi´on que representa el balance global del flujo del gradiente, donde Z φO = − ρ(r)dΩO (2.60) Ωo

y φd =

Z

Ωd

δ(r − rd )dΩd = 1

(2.61)

Con las funciones superarm´onicas hay m´as posibilidades de jugar con el balance del flujo del gradiente para conseguir un distribuci´on del mismo m´as homog´eneo. Las posibles condiciones de contorno son las de Dirichlet (u = g en ∂Ω), ∂u ∂u = g) y las mixtas (αu + β ∂n = g). El tratamiento las de Neuman ( ∂n detallado con la demostraci´on de soluci´on u ´nica se muestra en los apartados del ap´endice A.3.3, A.3.4 y A.3.5. La posibilidad de combinar las condiciones de contorno de Dirichlet y de Neuman, de forma que haya intervalos del mismo donde se aplique una u otra, permite controlar la distribuci´on del flujo (derivada direccional) en distintos tramos del contorno y por consiguiente en el domino. Considerando el modelo:  2  ∇ u = δ (r − rd ) ∀r ∈ Ω 

∂u α (ξ) u + β (ξ) ∂n = h ∀ξ ∈ ∂Ω

donde ξ parametriza el contorno y h es una funci´on gen´erica. Realizando

´ CAP. 2. FUNCIONES ARMONICAS

31

Figura 2.11: Condiciones de contorno mixtas. una partici´on del contorno en funci´on de la condici´on fijada (figura 2.11) ∂Ω = ∂ΩD + ∂ΩN ,

(2.62)

la expresi´on del flujo neto en todo el contorno ser´a, teniendo en cuenta el balance global del flujo del gradiente: I I I I ∂u ∂u ∂u ∂u d∂Ω = d∂Ω + d∂Ω = d∂Ω + 0 = 1, (2.63) ∂Ω ∂n ∂ΩD ∂n ∂ΩN ∂n ∂ΩD ∂n Con lo que es posible, mediante la gesti´on adecuada de las condiciones de contorno, localizar el “flujo disponible”(y por tanto el gradiente) en zonas convenientes, quiz´as alejadas del punto destino (apartado 3.4)

2.4.

C´ alculo anal´ıtico y obtenci´ on num´ erica de soluciones

Con la finalidad de ilustrar las propiedades de las funciones arm´onicas relativas a la planificaci´on de movimientos en rob´otica, es adecuado realizar el c´alculo anal´ıtico de algunos modelos b´asicos. La utilizaci´on de las funciones de Green es id´onea para este fin, ya que permite discernir entre las distintas partes que afectan a la soluci´on y completarla mediante una simple suma. La utilizaci´on del teorema de Green permite aislar la influencia de las condiciones de contorno sobre la funci´on potencial y sus propiedades en el dominio del espacio de configuraciones, la proposici´on de distintos modelos de C-espacios con obst´aculos regulares y las condiciones de compatibilidad ser´an la base para la exploraci´on de propiedades de la funci´on potencial ante distintas situaciones de contorno. En el ap´endice A.4 se muestra, para consulta, el teorema de Green. Funci´ on de Green Se define la funci´on de Green, o respuesta impulsional, como la soluci´on de ∇2 G = δ(r − rp ),

(2.64)

´ CAP. 2. FUNCIONES ARMONICAS

32

Figura 2.12: Contorno rectangular y funci´on delta. con r = (q1 , q2 , q3 , . . . , qn ), condiciones de contorno homog´eneas (G = 0) y δ(r − rp ) funci´on delta situada en rp . La funci´on de Green ser´a funci´on de la variable posici´on r y del punto de situaci´on de la funci´on delta rp , G(r, rp ). En general, la posici´on de la funci´on delta ser´a considerada como variable, dependiendo del modelo a resolver.

2.4.1.

Modelo base. An´ alisis de la influencia de las condiciones de contorno mediante las Funciones de Green

El modelo b´asico (sin considerar obst´aculos) est´a constituido por un C-espacio con dominio rectangular donde se define una funci´on potencial u que cumple con la ecuaci´on de Poisson ∇2 u = δ(r − rd ), con la funci´on delta situada en el punto objetivo destino de la trayectoria. Para aislar en la soluci´on la influencia de las condiciones de contorno se utilizan las funciones de Green; de manera que en el mismo dominio se define, adem´as de la funci´on potencial u, otra funci´on G con la que se completa el modelo:

∇2 G = δ(r − rp ), con condiciones de contorno homog´ eneas ∇2 u = δ(r − rd ), con condiciones de contorno no homog´ eneas

  

Donde G(r, rp ) es una funci´on de Green gen´erica, soluci´on de ∇2 G = δ(r − rp ), donde rp ser´a un punto considerado como variable, en contraste con g(r−rd ), que es la funci´on soluci´on de ∇2 g = δ(r − rd ), donde rd es el punto fijo de posici´on de la funci´on delta que marcar´a el m´ınimo global objetivo. Utilizando la 2a ecuaci´on de Green (ecuaci´on A.76) con las funciones u y

´ CAP. 2. FUNCIONES ARMONICAS G(r, rp ), se obtiene: Z

33

 I  ∂G ∂u u∇ G − G∇ u dVp = u −G dSp ∂n ∂n S V  Z I  ∂G ∂u (uδ(r − rp ) − G(r, rp )δ(r − rd )) dVp = u −G dSp ∂n ∂n V S  I  ∂u ∂G −G dSp u(rp ) − G(rd , rp ) = u ∂n ∂n S 2

2




por otro lado, las funciones de Green cumplen con la propiedad G(rd , rp ) ≡ G(rp , rd ) y G(rp , rd ) ≡ g(rp , rd ), con lo cual I I ∂u ∂G − G dS, (2.65) u(rp ) = g(rp , rd ) + u ∂n S S ∂n donde, en el segundo miembro de la igualdad, el primer t´ermino representa la influencia de la respuesta impulsional, el segundo representa la influencia de las condiciones de contorno de Dirichlet y el tercero la influencia de las condiciones de contorno de Neuman. Sin embargo, puesto que G = 0 en el contorno ∂Ω, se tiene: I ∂G u(rp ) = g(rp , rd ) + u (2.66) dS S ∂n El hecho de que G valga cero en el contorno (condiciones de contorno homog´eneas) facilita y hace posible la obtenci´on de la funci´on de Green para los modelos sencillos estudiados, pero el t´ermino correspondiente a las condiciones de contorno de Neuman queda anulado. Por este motivo, su estudio se realizar´a mediante el uso de las series de Fourier en el apartado 3.4. Con la utilizaci´on de funciones potenciales superarm´onicas el modelo base presentado es  ∇2 G = δ(r − rp ), con condiciones de contorno homog´ eneas  ∇2 u = δ(r − rd ) − ρ(r), con condiciones de contorno no homog´ eneas

Y la expresi´on obtenida de la 2a ecuaci´on de Green es: Z I ∂G u(rp ) = g(rp , rd ) + GρdV + u dS V S ∂n Igual que para las funciones potenciales arm´onicas el t´ermino representa la influencia de las condiciones de contorno.

2.4.2.



(2.67) H

S

dS u ∂G ∂n

Obtenci´ on de soluciones mediante c´ alculo num´ erico

La obtenci´on de soluciones utilizando las ecuaciones 2.66 y 2.67 es factible mediante la utilizaci´on de los m´etodos num´ericos de integraci´on cl´asicos y un

´ CAP. 2. FUNCIONES ARMONICAS

34

programa de c´alculo como el MATLAB. El contorno y los obst´aculos deben estar constituidos por formas geom´etricas regulares globales o por tramos de forma que sea posible parametrizar y definir variables de integraci´on que recorran el contorno y el volumen para posteriormente realizar la integraci´on num´erica. Puesto que el objetivo es obtener soluciones que sirvan para representar gr´aficamente la influencia de las condiciones de contorno y de la funci´on ρ, se plantean modelos b´asicos sencillos en 2D. En este caso la soluci´on gen´erica es Z I ∂G u(rp ) = g(rp , rd ) + GρdS + u dl (2.68) S C ∂n Para poder obtener soluciones mediante procedimientos num´ericos es necesario parametrizar la superficie y la curva cerrada de integraci´on Z b I ∂G ′ ∂G dl = |r (t)|dt, (2.69) u u ∂n a C ∂n donde

∂G = f [r(t)] (2.70) ∂n Para poder calcular num´ericamente la integral de contorno, primero se discretiza el integrando como funci´on del par´ametro t y despu´es se calcula la integral utilizando el m´etodo num´erico correspondiente (algoritmo 2.1) u

for t=1:TM y(t) = f [~r(t)] end; Io = trapz(t, y) Algoritmo 2.1: C´alculo num´erico de la integral de contorno mediante MATLAB

2.5.

Aportaci´ on

En este cap´ıtulo se ha realizado un estudio taxon´omico de las propiedades de las funciones arm´onicas aplicadas a la generaci´on de movimientos en rob´otica, diferenciando aquellas que son favorables de las que plantean una seria limitaci´on para su uso como funciones potenciales. Asimismo el cap´ıtulo presenta la metodolog´ıa para el c´alculo anal´ıtico de las funciones arm´onicas con el objetivo de validar las aportaciones que se har´an de cara a su mejor adecuaci´on a la planificaci´on de movimientos. - Una funci´on potencial debe tener como principal cualidad, para la obtenci´on de trayectorias en el espacio de configuraciones, la ausencia de m´ınimos locales en el dominio en el que est´a definida. Esta cualidad determina que un planificador de trayectorias basado en funciones arm´onicas ser´a completo; es

´ CAP. 2. FUNCIONES ARMONICAS

35

decir siempre encontrar´a una soluci´on si ´esta existe. - Las trayectorias obtenidas mediante el uso de las funciones arm´onicas ser´an suaves: es decir, la curvatura de las curvil´ıneas obtenidas como proyecci´on sobre el espacio de configuraciones por el seguimiento del gradiente de la funci´on arm´onica, ser´an ´optimamente reducidas para cualquier tipo de entorno. - Para determinadas condiciones de contorno la ecuaci´on de Poisson planteada, sobre cualquier entorno posible del C-espacio, tendr´a soluci´on y ´esta ser´a u ´nica. De otra manera, para que exista una soluci´on, las condiciones de contorno planteadas (valor y derivada de la funci´on) no pueden ser cualesquiera; debiendo de cumplir estas unas determinadas condiciones denominadas condiciones de solubilidad. - Variando ligeramente las condiciones de contorno la soluci´on experimenta, tambi´en, una ligera variaci´on. Esta propiedad permite la aplicaci´on de las funciones arm´onicas en entornos cambiantes, as´ı como la implementaci´on de algoritmos en los que la exploraci´on del C-espacio se va realizando sucesivamente. - La distribuci´on del gradiente, en principio, ser´a desigual; obteni´endose zonas donde este tomar´a valores extremadamente bajos y obligando (si el punto inicial se encuentra en ellas) a tiempos de computaci´on inadmisiblemente altos, especialmente si el n´ umero de grados de libertad es elevado, lo que har´ıa inviable el m´etodo. Ser´a necesario, por tanto, desarrollar un m´etodo que contrarreste este problema.

36

´ CAP. 2. FUNCIONES ARMONICAS

Cap´ıtulo 3 Alternativas de moldeado de las funciones arm´ onicas para su uso en la planificaci´ on de movimientos Teniendo en cuenta las condiciones de solubilidad que debe de cumplir todo modelo, la aplicaci´on de condiciones de contorno no uniformes de Dirichlet y mixtas suponen un recurso de moldeado de la funci´on potencial como soluci´on de un modelo formado por la ecuaci´on de Poisson y sus condiciones de contorno. La utilizaci´on de subdominios en los cuales la funci´on potencial ser´a definida como arm´onica, superarm´onica y subarm´onica, es el otro recurso utilizable para la obtenci´on de una funci´on potencial con una distribuci´on de gradiente m´as eficiente. Con la finalidad de solucionar el problema de las regiones planas o bajo gradiente se propone la utilizaci´on de las funciones subarm´onicas y superarm´onicas, juntamente con el establecimiento adecuado de las condiciones de contorno. En este cap´ıtulo se cumplen los siguientes objetivos: 1. Estudiar el efecto de las condiciones de contorno de Dirichlet sobre el moldeado de la funci´on potencial (apartado 3.3). 2. Estudiar el efecto de las condiciones de contorno de Dirichlet no uniformes sobre el moldeado de la funci´on potencial (apartado 3.4). 3. Estudiar el efecto de las condiciones de contorno mixtas (DirichletNeuman) sobre el moldeado de la funci´on potencial (apartado 3.5). 4. Estudiar el efecto de la utilizaci´on de las funciones sub y superarm´onicas, por subregiones, en el moldeado de la funci´on potencial (apartado 3.6). 37

38

CAP. 3. ALTERNATIVAS DE MOLDEADO

3.1.

Modelo gen´ erico

En este punto, es conveniente plantear un modelo gen´erico que contemple todas las posibilidades de modelado de los obst´aculos y del contorno, con la finalidad de estudiar e ilustrar las propiedades en cada caso. El modelo gen´erico planteado se aplica sobre un dominio Ω particionado en tres subdominios o regiones: regi´on subarm´onica ΩP , regi´on superarm´onica ΩO y regi´on homog´enea ΩH ; con un contorno particionado en dos zonas: zona donde se aplica las condiciones de contorno de Dirichlet ∂ΩD y zona donde se aplica las codiciones de contorno de Neuman ∂ΩN  2  ∇ u = ρ(r) + δ (r − rd ) en Ω 

∂u αu + β ∂~ = h en ∂Ω n

donde δ (r − rd ) modela el punto objetivo situado en rd , ρ(r) modela las regiones sub y superarm´onicas:  −ρO ∀r ∈ ΩO y ρO > 0      ρP ∀r ∈ ΩP y ρP > 0 ρ(r) =      0 ∀r ∈ / (ΩO ∪ ΩP ) y

  α = 1, β = 0 ∀r ∈ ∂ΩD 

α = 0, β = 1 ∀r ∈ ΩN

Aplicando la 2a ecuaci´on del teorema de Green (tal como se hace en el apartado 2.4.1)  Z I  ∂G ∂u u u(rp ) = g(rp , rd ) + ρGdV + −G dS, (3.1) ∂~n ∂~n V S

donde:

G es la funci´on de Green o respuesta impulsional del modelo:  2  ∇ G = δ (rp − rd ) ∀rd ∈ Ω 

condiciones de contorno homog´ eneas en ∂Ω

g(r, rd ) es la respuesta a la funci´on delta situada en el punto objetivo: g(r, rd ) = G(r, rd )

(3.2)

La integral de volumen determina la influencia de las regiones sub y superarm´onicas.

CAP. 3. ALTERNATIVAS DE MOLDEADO

39

La integral de superficie representa las condiciones de contorno: Condiciones de contorno de Dirichlet:− Condiciones de contorno de Neuman:

H

H

S

S


u ∂G dS ∂~ n


∂u dS G ∂~ n

Como se aprecia en la ecuaci´on 3.1, la parte fundamental del m´etodo es la obtenci´on de la funci´on de Green, lo cual s´olo es posible para modelos sencillos. De todas formas el objetivo es ilustrar las propiedades de las funciones arm´onicas, sub y superarm´onicas, as´ı como la aplicaci´on de distintas condiciones de contorno, y esto es posible mediante el estudio de modelos b´asicos. Para poder analizar las distintas condiciones de contorno, la obtenci´on de la funci´on de Green se realizar´a mediante el m´etodo de las cargas imagen (condiciones de contorno de Dirichlet) y mediante las series de Fourier (condiciones de contorno de Dirichlet y de Neuman), en ambos casos se plantear´an modelos representativos sencillos y se obtendr´a resultados utilizando m´etodos num´ericos mediante el programa MATLAB.

3.2.

Condiciones de contorno de Dirichlet. Obtenci´ on de la funci´ on de Green mediante el m´ etodo de las cargas im´ agenes, aplicado a funciones de dos variables

Mediante el m´ etodo de las im´ agenes se puede obtener la funci´on de Green, en un dominio con condiciones de contorno homogeneas, utilizando la funci´on de Green de todo el espacio; las sucesivas cargas imagen son r´eplicas especulares ficticias, que se forman fuera del dominio al considerar el contorno como si fuera un espejo. De esta forma se puede explorar las propiedades de la soluci´on de forma m´as intuitiva y sacar conclusiones respecto de la relaci´on entre la distribuci´on del gradiente y las condiciones de contorno. Para que el m´etodo sea factible, y produzca soluciones observables, se aplica a dominios en dos dimensiones.

3.2.1.

Funci´ on de Green para todo el espacio

De la ecuaci´on 2.50 se obtiene la funci´on de Green para todo el espacio, para un espacio de n dimensiones:   ∂g 1 (r − rd ) 1 (r − rd ) = = (3.3) · ∂(r − rd ) S |r − rd | |r − rd | S g(r, rd ) =

Z

dR S

(3.4)

40

CAP. 3. ALTERNATIVAS DE MOLDEADO

donde R = |r − rd |. Para n = 2, se tiene: Z 1 dR = ln R + C g(r, rd ) = 2πR 2π

(3.5)

y para n > 2, teniendo en cuenta 2.51:
Z
Γ n2 + 1 Γ n2 + 1 dR g(r, rd ) = =− +C n n Rn−1 nπ 2 n(n − 2)π 2 Rn−2

(3.6)

para un espacio infinito, en ambos casos.

3.2.2.

Ejemplos basados en condiciones de contorno homog´ eneas: Rect´ angulo y circunferencias obst´ aculo

Dominio semiplano. Aplicando el m´etodo a la obtenci´on de la funci´on de Green en un dominio definido por el semiplano y < 0, con condiciones de contorno homog´eneas en y = 0, obtenemos la carga imagen tal como se muestra en la figura 3.1a. La funci´on buscada es la soluci´on de ∇2 G = δ 0 − δ I

(3.7)

aplicada a todo el espacio. De acuerdo con el principio de superposici´on, la respuesta es la suma de las dos respuestas individuales: G=

1 1 1 |r − r0 | ln |r − r0 | − ln |r − rI | = ln 2π 2π 2π |r − rI |

De manera que para y = 0, |r − rI | = |r − r0 | y G =

1 2π

(3.8)

ln 1 = 0

Dominio franja. Para un dominio determinado por la franja 0 < y < L, las cargas imagen de δ0 se obtienen considerando dos espejos situados en y = 0 y y = L, tal como se muestra en la figura 3.1b. De esta forma se obtienen infinitas im´agenes cuyo signo de carga es alternado consecutivamente y el modelo aplicado a todo el espacio es: ∇2 G =

∞ ∞ X X (−1)(i+1) δ2i+1 (−1)i δ2i + i=0

(3.9)

i=0

cuya soluci´on es: G=

∞ ∞ 1 X 1 X 1 = (−1)i ln |r2i | + (−1)i ln 2π i=0 2π i=0 |r2i+1 | ∞ |r2i | 1 X (3.10) (−1)i ln 2π i=0 |r2i+1 |

CAP. 3. ALTERNATIVAS DE MOLDEADO

41

Figura 3.1: M´etodo de las im´agenes para la obtenci´on de la funci´on de Green: a) Semiplano, b) Franja de manera que para y = L, se cumple que |r2i+1 | = |r2i | ∀i|i = 0, 1, 2, · · · por lo que

∞ 1 X (−1)i ln 1 = 0 G= 2π i=0

(3.11)

(3.12)

Tambi´en se puede organizar de otra manera: G=

∞ 1 |r0 | 1 X |r2i | ln + (−1)i ln 2π |r2 | 2π i=2 |r2i−3 |

(3.13)

de forma que para y = 0 tambi´en se cumple que |r2 | = |r0 | y

|r2i−3 | = |r2i | ∀i|i = 2, 3, 4, · · ·

(3.14)

y tambi´en se cumplen las condiciones de contorno homog´eneas: ∞ 1 1 X G= ln 1 + (−1)i ln 1 = 0 2π 2π i=2

(3.15)

La posici´on de las sucesivas cargas im´agenes presentan una ley de recurrencia tal como se aprecia en la figura 3.1b, de forma que es posible obtener una expresi´on como un sumatorio de una serie. ∞ [(x − x0 )2 + (y − 2iL − y0 )2 ] 1 X (−1)i ln G(x, y) = 4π i=0 [(x − x0 )2 + (y + 2iL + y0 )2 ]

(3.16)

42

CAP. 3. ALTERNATIVAS DE MOLDEADO

Figura 3.2: M´etodo de las im´agenes para la obtenci´on de la funci´on de Green. Dominio rectangular Los t´erminos correspondientes a las cargas im´agenes m´as alejadas tienen un influencia despreciable y se puede prescindir de ellos. Es decir, conforme i · L es m´as grande el t´ermino del sumatorio correspondiente es m´as peque˜ no, cumpli´endose que l´ım ln

iL→∞

[(x − x0 )2 + (y − 2iL − y0 )2 ] =0 [(x − x0 )2 + (y + 2iL + y0 )2 ]

(3.17)

Dominio rect´ angulo. Un dominio definido por un rect´agulo delimitado por   0 < x < Lx 

0 < y < Ly

se puede considerar como la intersecci´on de dos franjas, y las cargas imagen de δ0 se pueden obtener como las r´eplicas sobre los espejos situados en x = 0, x = Lx , y = 0 y y = Ly (figura 3.2). La posici´on de las sucesivas cargas im´agenes presentan una ley de recurrencia, de forma que es posible obtener una expresi´on como un doble sumatorio de una serie: 2

∇ G=

∞ ∞ X X

(−1)

j=0 i=0 2

∇ G=

(i+j)

δ(2i,2j) −

∞ ∞ X X j=0 i=0

∞ ∞ X X

(−1)(i+j) δ(2i+1,2j+1)

(3.18)

j=0 i=0

  (−1)(i+j) δ(2i,2j) − δ(2i+1,2j+1)

(3.19)

CAP. 3. ALTERNATIVAS DE MOLDEADO

43

Figura 3.3: Funci´on potencial de Green con cargas im´agenes para un dominio rectangular: a) Funci´on de Green, b) Funci´on potencial de las cargas imagen, c) Situaci´on de las cargas imagen. cuya soluci´on es: ∞ ∞ |r − r(2i,2j) | 1 XX G= (−1)(i+j) ln 2π j=0 i=0 |r − r(2i+1,2j+1) |

(3.20)

de manera que para x = 0, x = Lx , y = 0, y = Ly , se cumple que |r − r(2i+1,2j+1) | = |r − r(2i,2j) | ∀i, j|i, j = 0, 1, 2, · · ·

(3.21)

por lo que en el contorno: ∞ ∞ 1 XX G= (−1)(i+j) ln 1 = 0 2π j=0 i=0

(3.22)

Los t´erminos correspondientes a las cargas im´agenes mas alejadas tienen un influencia despreciable, siendo nula en el infinito; por lo que se puede tener una aproximaci´on razonable manteniendo unos pocos t´erminos. La ecuaci´on 3.20 puede expresarse en t´erminos de x e y como: ∞ ∞ 1 XX [(x − 2jLx − x0 )2 + (y − 2 · i · Ly − y0 )2 ] (i+j) G(x, y) = (−1) ln 4π j=0 i=0 [(x − 2 · j · Lx − x0 )2 + (y + 2iLy + y0 )2 ]

(3.23) Utilizando la ecuaci´on 3.1, para ρ = 0 y G dada por la ecuaci´ın 3.23, mediante el MATLAB y sus funciones num´ericas, se obtiene la soluci´on que se representa en la figura 3.3.

Dominio exterior a una circunferencia. La obtenci´on de la funci´on de Green por el m´etodo de las im´agenes se realiza de acuerdo con el modelo de la figura 3.4. Esta funci´on ser´a la soluci´on de ∇2 G = δ 0 − δ I

(3.24)

Considerando la circunferencia centrada en el origen y de acuerdo con el principio de superposici´on, la respuesta es la suma de las dos respuestas

44

CAP. 3. ALTERNATIVAS DE MOLDEADO

Figura 3.4: Contorno infinito, funci´on delta y obst´aculo circular. Carga imagen. individuales: 1 1 |r − r0 |2 1 ln |r − r0 | − ln |r − rI | + C = ln + C (3.25) G= 2π 2π 4π |r − rI |2 De manera que para |r| = R, G debe valer cero. Lo cual ser´a cierto si se cumple que |r − r0 |2 = k|r − rI |2 , (3.26)

sobre la circunferencia y C = − (1/4π) ln k. Se demuestra que esto se cumple para una posici´on de la imagen sobre la l´ınea radial que une el centro de la circunferencia y la carga, tal que: rI =

r0 R2

(3.27)

k=

R2 r02

(3.28)

y

As´ı, la funci´on de Green ser´a:   |r − r0 |2 R2 1 ln G= 4π |r − rI |2 r02 y aplicando el teorema de los cosenos, se tiene:  2 2  1 R r + r02 − 2rr0 cos ϕ G= ln 2 2 4π r0 r + rI2 − 2rrI cos ϕ Puesto que rI =

(3.29)

(3.30)

R2 , r0

 2  1 r2 + r02 − 2rr0 cos ϕ R G= ln 2 2 , 4π r0 r + R4 /r02 − 2r(R2 /r0 ) cos ϕ

(3.31)

donde ϕ es el ´angulo entre r y r0 y, donde se observa que para r = R se cumple que G = 0 (en la figura 3.5 se muestra la representaci´on gr´afica de esta funci´on).

CAP. 3. ALTERNATIVAS DE MOLDEADO

45

Figura 3.5: Funci´on potencial de Green para un dominio exterior a una circunferencia: a) Funci´on de Green, b) Funci´on potencial de las cargas, c) Situaci´on de las carga imagen.

Figura 3.6: Cargas imagen para la obtenci´on de la funci´on de Green en un dominio comprendido entre un rectangulo y una circunferencia interior. Dominio rect´ angulo con una circunferencia interior. En este modelo el contorno esta formado por un contorno exterior rectangular (dominio del espacio de configuraciones) y un contorno interior circular que representa un obst´aculo. Seg´ un la figura 3.6, las cargas imagen se forman por las distintas r´eplicas especulares. Excepto la primera carga imagen, en el interior de la circunferencia obst´aculo, todas las r´eplicas se concentran en un reducido espacio y, por ser el signo de estas alternado, se considera su efecto neutro o nulo. Por lo tanto, el modelo v´alido aproximado es 2

∇ G≃

∞ ∞ X X j=0 i=0

  (−1)(i+j) δ0(2i,2j) − δ0(2i+1,2j+1) + δI(2i+1,2j+1) − δI(2i,2j)

cuya soluci´on es, aplicando superposici´on:   ∞ ∞ |r0(2i,2j) ||rI(2i+1,2j+1) | 1 XX (i+j) (−1) ln G≃ 2π j=0 i=0 |r0(2i+1,2j+1) ||rI(2i,2j) |

(3.32)

(3.33)

46

CAP. 3. ALTERNATIVAS DE MOLDEADO

Figura 3.7: Funci´on potencial de Green para un dominio comprendido entre un rectangulo y una circunferencia interior: a) Funci´on de Green (incluyendo el potencial dentro del obst´aculo), b) Funci´on potencial de las cargas, c) Situaci´on de las cargas imagen. de manera que para x = 0, x = Lx , y = 0, y = Ly , se cumple que |r0(2i+1,2j+1) | = |r0(2i,2j) | y

|rI(2i+1,2j+1) | = |rI(2i,2j) |;

∀i, j|i, j = 0, 1, 2, · · · (3.34)

por lo que en el contorno: ∞ ∞ 1 XX (−1)(i+j) ln 1 = 0 G≃ 2π j=0 i=0

(3.35)

Por su parte tambi´en G ≃ 0 en la circunferencia; ya que la influencia de las cargas especulares exteriores al rect´angulo inducir´an, a su vez, im´agenes en el interior a la circunferencia aglutinadas cerca del centro, cancel´andose por tanto unas a otras. La funci´on 3.33 puede expresarse en t´erminos de x e y como:  ∞ ∞ 1 XX [(x − 2jLx − x00 )2 + (y − 2iLy − y00 )2 ] (i+j) G(x, y) ≃ (−1) ln · 4π j=0 i=0 [(x − 2jLx − x00 )2 + (y + 2iLy + y00 )2 ] [(x − 2jLx − x0I )2 + (y + 2iLy + y0I )2 ] [(x − 2jLx − x0I )2 + (y − 2iLy − y0I )2 ]



(3.36)

La programaci´on, mediante MATLAB, de la soluci´on se representa en la figura 3.7, en ella se aprecia el cumplimiento de las suposiciones y del resultado, tanto de las cargas imagen, de las condiciones de contorno homog´eneas, como de la funci´on de Green obtenidas.

3.2.3.

Contorno infinito, funci´ on delta y obst´ aculo circular

En este caso se puede considerar el obst´aculo como el contorno. Es decir, existir´a un contorno interno circular y la integral de circulaci´on del contorno

CAP. 3. ALTERNATIVAS DE MOLDEADO

47

Figura 3.8: Obst´aculo circular. ser´a la de la circunferencia interna (figura 3.8), ya que el gradiente de G en el infinito es cero. De esta forma, la expresi´on global del potencial ser´a, para condiciones de contorno gen´ericas de Dircichlet: I ∂G dS (3.37) u=g+ f (r) ∂n O Para un espacio de 2 dimensiones: Dada la ecuaci´on 3.31:  2  R 1 r2 + r02 − 2rr0 cos ϕ ln 2 2 G= , 4π r0 r + R4 /r02 − 2r(R2 /r0 ) cos ϕ

(3.38)

para evaluar la integral de contorno, que determina la influencia de las condiciones de contorno, se obtiene la derivada de la funci´on de Green con respecto a r0 en el contorno: R ∂G 1 − (r/R)2 = (3.39) ∂r0 r0 =R 2π r2 + R2 − 2Rr cos ϕ

Por lo tanto la aportaci´on de las condiciones de contorno (f (r0 )) a la soluci´on de la ecuaci´on en derivadas parciales viene dada por la integral: Z I 1 − (r/R)2 ∂G R 2π f (θ ) Rdθ0 IO (r, θ) = f (r0 ) dL = 0 ∂r0 r0 =R 2π 0 r2 + R2 − 2Rr cos(θ − θ0 ) O (3.40) Z 2π 1 R2 − r 2 IO (r, θ) = f (θ0 ) 2 dθ0 (3.41) 2π 0 r + R2 − 2Rr cos(θ − θ0 ) denominada integral de Poisson. A una distancia elevada del contorno circular, la influencia de sus condiciones de contorno se uniformiza a un valor

48

CAP. 3. ALTERNATIVAS DE MOLDEADO

Figura 3.9: Condiciones de contorno constante (uC = KC = 5) para un obst´aculo circular. constante: 1 l´ım IO (r, θ) = − r→∞ 2π

Z

2π

f (θ0 )dθ0 = −fM ED (θ0 ),

0

(3.42)

que es el valor medio de las condiciones de contorno y que es independiente de la posici´on del punto como funci´on de θ. Aplicando el teorema del valor medio para integrales del producto de dos funciones, para n dimensiones en general se cumple: H H ∂G
H S f (rC )dS dS = f (r ) dS = u · 1 = IOM ED = S f (rC ) ∂G M ED C medio ∂n S S ∂n IOM ED =

H

S

f (rC )dS S

(3.43)

Adem´as, teniendo en cuenta que u(r, θ) = g(r, θ) + Io(r, θ), se cumple: I ∂Io dS = 0 (3.44) O ∂r

para cualquier superficie cerrada que envuelva el obst´aculo, lo que quiere decir que: • El balance de la derivada de la funci´on influencia Io, se compensar´a haciendose cero. Es decir, si la pendiente de salida del obst´aculo tiene un valor elevado en alguna zona, puede provocar que la pendiente entrante en otra determine m´ınimos locales sobre el contorno. Condiciones de contorno: valor constante de la funci´ on potencial sobre el contorno de un obst´ aculo circular Para uC = KC , la influencia que tiene sobre la soluci´on viene dada por I Z ∂G 1 − (r/R)2 KC R 2π IO (r, θ) = KC Rdθ0 = KC dL = 2π 0 r2 + R2 − 2Rr cos(θ − θ0 ) O ∂r0 r0 =R (3.45)

CAP. 3. ALTERNATIVAS DE MOLDEADO

49

resultado coherente con las condiciones de solubilidad, ya que I ∂G KC dL = KC · 1 = KC O ∂r0 r0 =R

(3.46)

El valor constante de la integral determina que:

• La influencia, de un valor contante de las condiciones de contorno, sobre el gradiente de la funci´on potencial, es nula en todo el dominio y solo afecta al valor de la funci´on en una constante sumatoria, figura 3.9.

3.3.

Moldeado de la funci´ on potencial mediante la aplicaci´ on de condiciones de contorno no uniformes de Dirichlet

La aplicaci´on de condiciones de contorno de Dirichlet no uniformes tiene un efecto de moldeado sobre la funci´on potencial soluci´on que es necesario conocer. El contorno gen´erico a tener en cuenta, est´a formado por un contorno rectangular externo ∂ΩR y un contorno interno definido por los obst´aculos ∂ΩO . Sobre el dominio limitado por este contorno se aplica la ecuaci´on de Poisson ∇2 u = δ(rd ). Las condiciones de contorno no uniformes de Dirichlet suponen que el valor de la funci´on potencial en el contorno ser´a funci´on del punto considerado:  u(r)|C = f (rC )  H

∂u dS n S ∂~

=1



Aplicando la segunda ecuaci´on de Green se obtiene la soluci´on: I ∂G dS u = g + f (r) ∂n S

(3.47)

Inicialmente se estudian varios modelos b´asicos y se extraen las conclusiones generales que se tendr´an en consideraci´on a la hora de dise˜ nar el algoritmo planificador.

3.3.1.

Dominio en forma de franja

La funci´on de Green viene dada por la expresi´on 3.16, donde el t´ermino gen´erico del sumatorio es: Gi (x, y; xp , yp ) =

1 (x − xp )2 + (y − 2iL − yp )2 ln 4π (x − xp )2 + (y + 2iL + yp )2

(3.48)

50

CAP. 3. ALTERNATIVAS DE MOLDEADO

y sus derivadas en la direcci´on normal del contorno, yp = L y yp = 0: ∂Gi (x, y; xp , yp ) = ∂yp yp =L

  2(y − (2i + 1) L) 2(y + (2i + 1) L)) 1 , + − 4π (x − xp )2 + (y − (2i + 1) L)2 (x − xp )2 + (y + (2i + 1) L)2 (3.49)

∂Gi (x, y; xp , yp ) = ∂(−yp ) yp =0

  2(y − 2 · i · L) 2(y + 2 · i · L) 1 (3.50) + 4π (x − xp )2 + (y − 2 · i · L)2 (x − xp )2 + (y + 2 · i · L)2

De esta forma se obtiene la integral de contorno: Z +∞ I Z +∞ ∞ ∞ X X ∂Gi ∂G ∂Gi i i f (rp ) dxp + (−1) dxp f (r) f (rp ) (−1) dS = ∂n ∂yp yp =L ∂(−yp ) yp =0 −∞ S −∞ i=1 i=1

(3.51)

Condiciones de contorno dependientes de la distancia, para un dominio en franja: uC |y=L = uC |y=0 = k · |x − x0 | En este modelo el valor de la funci´on potencial en el contorno es proporcional a la distancia respecto de la situaci´on de la funci´on delta. I ∂G f (r) dS = ∂n S k

∞ X i=1

(−1)

i

Z

+∞ −∞

Z +∞ ∞ X ∂Gi ∂Gi i dxp +k dxp (−1) |x−x0 | |x−x0 | ∂yp yp =L ∂(−y ) p −∞ y =0 p i=1

(3.52)

En la figura 3.10 se muestran el modelo y los resultados. Se observa que es posible modular la derivada tangencial en el contorno, mientras que el balance global de la derivada normal en el mismo permanece constante. En este caso la derivada normal aportada por las condiciones de contorno es cero en todo el contorno y por lo tanto no existe un flujo de paso interferente, sin embargo: • Se establece una derivada tangencial que garantiza la existencia de un gradiente elevado incluso en los puntos alejados del destino, lo cual es un efecto deseable.

CAP. 3. ALTERNATIVAS DE MOLDEADO

51

Figura 3.10: Condiciones de contorno dependiente de la distancia para un dominio franja.

Figura 3.11: Condiciones de contorno funci´on de la distancia sobre el per´ımetro de un obst´aculo circular: (a) Contorno del obst´aculo circular. (b) T´ermino correspondiente a la integral de contorno. (c) Funci´on potencial, se˜ nalando el m´ınimo sobre el contorno.

52

3.3.2.

CAP. 3. ALTERNATIVAS DE MOLDEADO

Condiciones de contorno proporcional con la distancia sobre el contorno de un obst´ aculo circular

Con la finalidad de establecer una derivada tangencial sobre el contorno obst´aculo, que determine un gradiente favorable en el entorno pr´oximo, se fija uC = k · l(θ0 ) de forma que la referencia se encuentre sobre una l´ınea que una el punto destino con el centro de la circunferencia y que intersecta la misma en los puntos P y Q (figura 3.11). La dependencia lineal tendr´a su valor m´aximo en el punto P y valdr´a cero en el punto Q I Z 2π k ∂G R2 − r 2 dL = IO (r, θ) = k·l(θ0 ) l(θ ) dθ0 0 2 ∂r0 r0 =R 2π 0 r + R2 − 2Rr cos(θ − θ0 ) O (3.53) • La influencia que tiene estas condiciones de contorno, sobre la distribuci´on del gradiente, es favorable ya que en conjunto se observa una cierta modulaci´on en la direcci´on del punto destino, sin embargo aparece la posibilidad (si se acent´ ua demasiado el afecto) de un m´ınimo local en el punto Q sobre el contorno.

3.4.

Moldeado de la funci´ on potencial mediante el uso de condiciones de contorno mixtas

En el caso de un dominio rectangular, con obst´aculos formando parte del contorno, es preferible utilizar las series de Fourier para obtener la funci´on de Green ya que nos permitir´a aplicar condiciones de contorno tanto de Dirichlet como de Neuman.

3.4.1.

Obtenci´ on de la funci´ on de Green utilizando las series de Fourier

Con la finalidad de poder estudiar la influencia que tienen las condiciones de contorno de Neuman en la distribuci´on del gradiente, se propone obtener la funci´on de Green soluci´on del siguiente modelo:  2 ∇ G = δ(x − x0 , y − y0 ) en Ω      G(0, y) = 0 y G(L, y) = 0, 0 0      ρP ∀r ∈ ΩP y ρP > 0 ρ(r) =      0 ∀r ∈ / (ΩO ∪ ΩP )

tal que ρO y ρP son valores constantes y positivos. Considerando la influencia, sobre la soluci´on, del obst´aculo Oi ∇2 ui = ρ(r),

(3.72)

tal que ρ(r) = y del pozo de potencial P

  −ρOi 

∀r ∈ ΩOi y ρOi > 0

0 ∀r ∈ / ΩOi ∇2 up = ρ(r),

(3.73)

tal que ρ(r) =

  ρP 

∀r ∈ ΩP y ρP > 0

0 ∀r ∈ / ΩP

y aplicando el principio de superposici´on, se obtiene la soluci´on influencia de los obst´aculos en un punto cualquiera del dominio u=

n X

ui + up

(3.74)

i=1

Similarmente, el gradiente en un punto ser´a la suma vectorial ∇u =

n X i=1

∇ui + ∇up

(3.75)

CAP. 3. ALTERNATIVAS DE MOLDEADO

59

Por otro lado, es conveniente homogeneizar la influencia de los obst´aculos en la soluci´on y en el gradiente resultante, de manera que la repulsi´on producida por un obst´aculo no dependa de su tama˜ no, para ello se define el concepto de masa de un obst´aculo como MO i = ρ O i V O i

(3.76)

de forma que, para un obst´aculo aislado, el flujo del gradiente que atraviesa una esfera que lo contiene es I Z φ Oi = ∇ui · dS = ∇2 udV = −ρOi VOi = −MOi (3.77) S

V

Para un obst´aculo aislado circular (2D), el valor m´aximo situado en el centro ser´ıa ρk 2 MO i uM AXOi = uCOi + RO = uCOi + , (3.78) i 4 4π Debi´endose cumplir, para uniformizar la influencia de los obst´aculos, que M O i = MO j = K M

∀(i, j) ∈ [n, n] y

| i 6= j

(3.79)

Lo que supone que el valor constante dado a la funci´on densidad (ρOi ) deber´a ser ponderado: KM (3.80) ρOi = V Oi El establecimiento de la regi´on parche en torno del objetivo debe de realizarse teniendo en cuenta la funci´on que debe desempe˜ nar: - Acentuar el gradiente en el dominio sin crear m´ınimos locales. Para cumplir este fin, la regi´on parche debe estar confinada por los obst´aculos, debe tener un tama˜ no reducido y su valor ρ no debe ser muy elevado. De este modo se minimiza el riesgo de m´ınimo local pero no se anula, con la modelizaci´on mediante coronas propuesta en la secci´on siguiente se elimina este riesgo. En el caso ideal del ejemplo en 2D, el parche es un c´ırculo centrado en el objetivo y los obst´aculos son tambi´en modelados como c´ırculos. Por u ´ltimo, el balance global (obst´aculos, parche y funci´on delta) del flujo, ser´a Z I ∇2 udV = 1 + MP − MO (3.81) ∇u · dS = φG = SG

VG

donde MP es la masa del parche, con signo contrario a la masa de los obst´aculos y valor MP = ρP Vp . Este flujo global φG , interesa que tenga un cierto valor positivo que deber´a aportar el contorno (en este caso situado en el infinito por simplicidad), de esta forma existir´a un gradiente m´ınimo en todo el dominio. Con la finalidad de que exista un determinado gradiente en el dominio que

60

CAP. 3. ALTERNATIVAS DE MOLDEADO

conduzca hacia el punto objetivo, el contorno debe aportar un gradiente definido por una Densidad Media de Flujo en el Contorno DM F C =

φG LCON T

=K

(3.82)

donde φG es el flujo global del gradiente aportado por el contorno, que cumple φG = 1 + MCP − MCO

(3.83)

Entonces la fijaci´on de φG condiciona el valor de ρP ρP =

φ G − 1 + MO MP = VP VP

(3.84)

En resumen y a la vista de los resultados, se concluye: • La definici´on de una zona ΩP , en torno de la funci´on delta, donde la funci´on ser´a subarm´onica, supone la generaci´on de un pozo de potencial que acent´ ua el gradiente en todo el dominio. • Manejando adecuadamente MP y MO se puede obtener un d´eficit positivo del balance del flujo entre los obst´aculos, el pozo de potencial y la funci´on δ; lo que determina la necesidad de aportaci´on de flujo del gradiente por parte del contorno y, por consiguiente, una mejora del mismo en el dominio. • El hecho de poder definir la masa de cada obst´aculo, hace el procedimiento m´as flexible para su combinaci´on con los m´etodos de exploraci´on aleatoria del espacio de configuraciones; en los cuales se obtiene informaci´on parcial del entorno, en etapas sucesivas (cap´ıtulo 6)

3.5.3.

Modelizaci´ on de los obst´ aculos mediante coronas circulares.

La mayor facilidad para determinar la superficie de contorno en lugar del volumen, tanto de los obst´aculos como del pozo de potencial, as´ı como una mejor definici´on de los contornos, sugiere la utilizaci´on de superficies de contorno en lugar de vol´ umenes para modelar los obst´aculos y el pozo de potencial. As´ı, utilizando la misma idea de modelizaci´on del apartado anterior, se realiza la modificaci´on de las regiones de vol´ umenes a coronas   ΩOi → ΩCOi →



ΩP

ΩCP



MCP = ρCP SCP

La definici´on de masa de las coronas ser´a, similarmente   MCOi = ρCOi SCOi

CAP. 3. ALTERNATIVAS DE MOLDEADO

61

Figura 3.16: Modelizaci´on de los obst´aculos y pozo de potencial mediante coronas circulares.

La modelizaci´on de un obst´aculo mediante una funci´on superarm´onica, en forma de corona circular, vendr´a dada por ∇2 u = ρ(r − rO ),

(3.85)

donde ρ tomar´a un valor constante negativo dentro de la corona circular, centrada en rO , y valdr´a cero en el resto   −ρk , ∀r ∈ R1 ≤| r − rO |≤ R2 ρ(r − rO ) =  0, ∀r ∈ (| r − rO |< R1 ) ∪ (| r − rO |> R2 ) Para r < R1

∇u = 0

Para R1 ≤ r ≤ R2

→

u = Constante

Sr ρk ∇u = −ρk =− Lr 2



R2 r− 1 r



(3.86)

(3.87)

y u = uR 2 −

Z

r ′

R2

∇u dr = uR2

 
ρk R12 r ρk 2 2 R2 − r − ln + 4 2 R2

(3.88)

El valor m´aximo del potencial se obtiene para r = R1    
ρk R12 ρk MCOi ρk R12 R1 R1 2 2 uM AXCOi = uR2 + ln = uR 2 + − ln , R2 − R1 − 4 2 R2 4π 2 R2 (3.89)

62

CAP. 3. ALTERNATIVAS DE MOLDEADO

donde la masa de la corona vale MCOi = ρCOi SCOi ≈ ρCOi eLCOi ,

(3.90)

ya que el espesor de la corona e ser´a bastante inferior al radio interior. De todas formas, como el objetivo es dar una cierta uniformidad a la influencia de los obst´aculos, con una aproximaci´on es suficiente. Para r > R2 ∇u = −

(R22 − R12 ) ρk 1 MCOi 1 · =− · 2 r 2π r

(3.91)

con el valor m´aximo del m´odulo del gradiente en r = R2 MCOi (R22 − R12 ) ρk =− (3.92) 2 2R2 2πR22 Imponiendo la condici´on de que la influencia de todos los obst´aculos sea la misma: MOi = KM , el valor ponderado de la funci´on ρ para la corona representativa del obst´aculo Oi ser´a ∇u = −

ρCOi =

KM /e LCOi

(3.93)

El producto escalar entre el gradiente y el vector normal, en cualquier punto del contorno, debe ser positivo (gradiente saliente del contorno) (∇u · n)|P > 0 ∀P ∈ Cont,

(3.94)

pero, aplicando condiciones de Dirichlet, el gradiente en un punto del contorno es la suma vectorial de los gradientes aportados por los obst´aculos ∇uCOi , por la corona pozo ∇uCP y por la funci´on delta ∇ud , de acuerdo con el principio de superposici´on. Para evitar que ∇u · n < 0 en uno o varios puntos del contorno, debe de cumplirse que DM F C >> 0. Pero esta condici´on no es suficiente, ya que puede ocurrir que un determinado obst´aculo est´e muy pr´oximo al contorno d(Oi , CON T )

<

λM IN

(3.95)

y su gradiente (saliente del mismo) supere el gradiente medio global del contorno, en este caso se producir´ıan m´ınimos locales en el contorno. Esta situaci´on se puede evitar de alguna de las maneras siguientes: Anexionando el obst´aculo al contorno.- Tiene la pega de que el pasillo entre el obst´aculo y el contorno queda eliminado. Aplicando al obst´aculo condiciones de contorno de Neuman.- Es siempre un buen recurso para aplicar en los casos problem´aticos y dificiles. Adem´as, si aplicamos Neuman a un obst´aculo, el balance del gradiente que tiene que aportar el contorno aumenta. Tiene la pega que las trayectorias pueden pasar rozando el obst´aculo.

CAP. 3. ALTERNATIVAS DE MOLDEADO

63

Figura 3.17: Adaptaci´on de la capa de pozo de potencial a un obst´aculo pr´oximo. Reducir ρ en el obst´aculo implicado y aumentarlo en la Corona del Pozo de potencial CP , para hacer que |∇ui | < |∇uM ED | Por lo tanto, en la utilizaci´on de las funciones sub y superarm´onicas para la planificaci´on de movimientos en rob´otica, la modelizaci´on mediante coronas o capas superficiales es id´onea porque: Permite mejorar la distribuci´on del gradiente en el dominio sin perder la definici´on de los contornos. La homogeneizaci´on de la influencia de los obst´aculos es m´as f´acil mediante la estimaci´on aproximada de la longitud de sus contornos. La utilizaci´on de una corona de atracci´on de potencial, en torno de la funci´on delta destino, mejora el gradiente sin el peligro de aparici´on de un m´ınimo local ya que, el gradiente aportado por la corona en su interior es nulo y su influencia nula en esa regi´on. As´ı, el modelo b´asicamente considera el contorno (al cual se le aplican condiciones de contorno homog´eneas de Dirichlet) y los obst´aculos que quedan definidos por regiones en forma de capas o coronas donde la funci´on potencial ser´a superarm´onica, adem´as una capa en torno del punto destino donde la funci´on ser´a subarm´onica que reforzar´a el gradiente de atracci´on hacia el objetivo. La obtenci´on de la regi´on capa denominada pozo de atracci´on CP , debe realizarse teniendo en cuenta el espacio libre disponible en torno del punto objetivo; estableciendo una distancia m´ınima a los obst´aculos dM IN y una cierta distancia m´axima al punto objetivo dM AX . As´ı, la capa pozo estar´a comprendida dentro de la corona, centrada en el punto objetivo, de radio

64

CAP. 3. ALTERNATIVAS DE MOLDEADO

m´ınimo dM IN y de radio m´aximo dM AX , sin que intersecte con los obst´aculos en ning´ un punto, tal como se muestra en la figura 3.17. Este m´etodo de representaci´on de obst´aculos, mediante funciones superarm´onicas, es presentado por I˜ niguez y Rosell en la publicaci´on [6] del apartado 7.1.

3.6.

Aportaci´ on

En este cap´ıtulo se ha realizado el estudio, con obtenci´on num´erica de resultados, de la influencia de las condiciones de contorno en las propiedades de la soluci´on de la ecuaci´on de Laplace. Igualmente se han analizado las propiedades de las funciones superarm´onicas y subarm´onicas. De esta manera, las aportaciones de este cap´ıtulo se resumen el los siguientes ´ıtems: - Se han desarrollado implementaciones para el c´alculo anal´ıtico de funciones arm´onicas en entornos simples en 2D: Utilizando las funciones de Green para resolver la ecuaci´on de Poisson en distintos modelos b´asicos del espacio de configuraciones, y utilizando el MATLAB como herramienta de c´alculo num´erico, se han calculado e ilustrado gr´aficamente distintas situaciones sobre las que se han podido extraer conclusiones fundamentales para el moldeado del gradiente. - Se ha propuesto el uso de condiciones de contorno de Dirichlet no uniformes como una opci´on para el moldeado del gradiente en las funciones arm´onicas. Bas´andose en esta posibilidad se ha desarrollado un m´etodo para la obtenci´on de funciones potenciales aplicables a la planificaci´on de movimientos; la posibilidad de aparici´on de un m´ınimo puntual en el contorno, as´ı como un determinado coste computacional, son los puntos d´ebiles del procedimiento. - Se ha propuesto la utilizaci´on de condiciones de contorno mixtas, con zonas donde se aplican condiciones de contorno de Dirichlet y zonas donde se aplican condiciones de contorno de Neuman, que permiten modelar el gradiente de la funci´on potencial resultante, de forma que se puede garantizar un gradiente m´ınimo en zonas alejadas del objetivo. - Se ha propuesto una opci´on novedosa, independientemente de las condiciones de contorno, para el moldeado del gradiente: La utilizaci´on de funciones superarm´onicas y subarm´onicas. El modelado de obst´aculos mediante funciones superarm´onicas, as´ı como la definici´on de una regi´on subarmonica en torno del punto objetivo, permite moldear el gradiente por efecto de balanceado de masas positivas y negativas. Bas´andose en estas propiedades se ha desarrollado otro m´etodo para la obtenci´on de funciones potenciales.

Cap´ıtulo 4 Planificaci´ on de movimientos en entornos din´ amicos y cambiantes En un entorno din´amico, constituido por obst´aculos en movimiento, es posible obtener una funci´on potencial dependiente del tiempo que sirva para la planificaci´on de trayectorias: Modelando el espacio de configuraciones mediante la ecuaci´on del calor, de forma que los obst´aculos m´oviles se correspondan con una funci´on generadora de calor, su soluci´on da una funci´on potencial con la que es posible predecir la trayectoria de un obst´aculo y crear los movimientos reactivos que eviten colisiones, tanto si el robot est´a en movimiento como si est´a en reposo. En el apartado 4.1 se realiza el an´alisis de un modelo matem´atico b´asico y en el 4.2 se deduce un m´etodo de control que corrije la posici´on o la trayectoria del robot, evitando colisiones con un obst´aculo en movimiento. Por otro lado, un m´etodo que utilizase una exploraci´on aleatoria incremental en un entorno est´atico; en cada paso exploraria un poco m´as del espacio de configuraciones, obteni´endose ligeras variaciones del mismo. Utilizando una funci´on potencial sobre este entorno cambiante, se realizar´ıan actualizaciones del c´alculo de dicha funci´on potencial: Cada nuevo c´alculo parte de un valor inicial que es calculado en el paso anterior. Por la propiedad de estabilidad de la soluci´ on, el nuevo valor de la funci´on potencial diferir´a poco del anterior y el tiempo de computaci´on ser´a muy reducido. En el apartado 4.3 se realiza un an´alisis te´orico de los casos presentados.

4.1.

Estudio anal´ıtico para condiciones iniciales y funci´ on obst´ aculo

Un entorno din´amico, en el espacio de configuraciones, formado por un obst´aculo en movimiento y contorno en el infinito, puede ser modelado por la ecuaci´on del calor; con una funci´on como foco generador de calor en movimiento, en representaci´on de un obst´aculo m´ovil, y unas determinadas 65

66

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

Figura 4.1: Funci´on influencia de difusi´on.

Figura 4.2: a) Funci´on generadora de calor, modelo de un obst´aculo en movimiento. b) Soluci´on de la funci´on obst´aculo en movimiento.

condiciones iniciales. La obtenci´on de una expresi´on anal´ıtica para la funci´on potencial soluci´on de la Ecuaci´on en Derivadas Parciales (EDP) y la representaci´on gr´afica correspondiente mediante MATLAB, permite extraer conclusiones que sirvan de base para el dise˜ no de la estrategia de control de trayectorias en entornos din´amicos. La complejidad del planteamiento requiere la consideraci´on de un modelo para una variable, con q(x, t) como funci´on generadora que representa el obst´aculo en movimiento:

  

∂u ∂t

2

= k ∂∂xu2 + q(x, t), −∞ < x < ∞

c.i. u(x, 0) = f (x)

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

67

Aplicando la transformada de Fourier F respecto de la variable x: R∞   F (ω) = −∞ f (x)e−jωx dx 

f (x) =

1 2π

R∞

−∞

F (ω)ejωx dω

se obtiene la Ecuaci´on diferencial Ordinaria (EDO)  ∂  ∂t U (ω, t) + kω 2 U (ω, t) = Q(ω, t) 

c.i. U (ω, 0) = F (ω)

cuya soluci´on ser´a la suma de la soluci´on homog´enea m´as la soluci´on particular U (ω, t) = Uh + Up (4.1) donde Uh = c(ω)e−kω

2t

(4.2)

y la soluci´on particular ser´a Up = e

−kω 2 t

Z

t

2

Q(ω, τ )ekω τ dτ,

(4.3)

0

tal como se puede demostrar sustituyendo en la EDO y teniendo en cuenta el primer principio de c´alculo. Por lo tanto, la soluci´on completa ser´a Z t 2 −kω 2 t −kω 2 t Q(ω, τ )ekω τ dτ, (4.4) U (ω, t) = c(ω)e +e 0

donde c(ω) = U (ω, 0) = F (ω)

(4.5)

Sustituyendo U (ω, t) = F (ω)e

−kω 2 t

+e

−kω 2 t

Z

t

2

Q(ω, τ )ekω τ dτ

(4.6)

0

Finalmente, aplicando el teorema de convoluci´on y teniendo en cuenta r π − x2 −kω 2 t −1 (4.7) e F e 4kt −−→ kt se obtiene la soluci´on de la EDP 1 u(x, t) = 2π

Z

∞

f (x) −∞

r

1 π − (x−x)2 e 4kt dx + kt 2π

Z tZ 0

∞

q(x, τ ) −∞

r

(x−x)2 π e− 4kt dxdτ k(t − τ ) (4.8)

68

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

Figura 4.3: Funci´on potencial y funci´on obst´aculo en movimiento. Tomando como par´ametro el tiempo o la distancia en la soluci´on mostrada en 4.2b. (a = 5m y v = 0,5m/s)

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

69

El primer t´ermino integral representa la influencia de las condiciones iniciales: r Z ∞ 1 π − (x−x)2 ud = (4.9) f (x) e 4kt dx 2π −∞ kt que se difunden dependiendo de k mediante la funci´on influencia: r π − (x−x)2 e 4kt g(x, x, t) = kt

(4.10)

representada en la figura 4.1. El segundo t´ermino integral representa la influencia de la funci´on q(x, t) generadora de calor: Z tZ ∞ r (x−x)2 π 1 (4.11) e− 4kt dxdτ q(x, τ ) uc = 2π 0 −∞ k(t − τ ) Suponiendo un obst´aculo en movimiento representado por la funci´on:   1, para x0 + vt − a2 < x < x0 + vt + a2 q(x, t) =  0, para el resto

donde x0 es la posici´on inicial del centro del obst´aculo, v su velocidad de desplazamiento y a su anchura (figura 4.2a) Entonces, el t´ermino integral (4.11) valdr´a Z t Z x0 +vτ + a r 2 (x−x)2 1 π uc = (4.12) e− 4kt dxdτ 2π 0 x0 +vτ − a2 k(t − τ ) Su representaci´on gr´afica, obtenida por m´etodos num´ericos con MATLAB, se muestra en las figuras 4.2b y 4.3 para a = 5u.d.l. y v = 0,5m/s.

4.2.

Control de la configuraci´ on y de la trayectoria de un robot en un entorno con obst´ aculos m´ oviles.

De acuerdo con los resultados cualitativos obtenidos en el apartado anterior, un obst´aculo en movimiento (representado por un foco de calor, tambi´en en movimiento) producir´a una funci´on potencial (soluci´on de la ecuaci´on del calor) que presentar´a una asimetr´ıa en la distribuci´on de las curvas de nivel (y en el gradiente) con respecto a dicho obst´aculo en movimiento: estas se concentrar´an en el frente del movimiento y se separar´an en la estela del mismo. El navegador que obtenga trayectorias a partir de esta funci´on potencial din´amica, deber´a evitar las colisiones mediante la deduci´on de movimientos correctores reactivos en respuesta al acercamiento de un determinado obst´aculo; para ello seguir´a la evoluci´on puntual de la funci´on potencial correspondiente a la posici´on actual del robot:

70

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

Figura 4.4: Funci´on potencial y obst´aculo en movimiento: v0 , velocidad del obst´aculo O. T1 , trayectoria prevista inicialmente. VT , velocidad prevista en el punto P . c, velocidad de acercamiento detectada en el punto P . vr , velocidad reactiva. Si un obst´aculo se acerca a un punto el valor de la funci´on aumentar´a en dicho punto y la variaci´on temporal relativa por cada coordenada (derivada temporal del gradiente), determinar´a el vector direcci´on de acercamiento o alejamiento (figura 4.4). Este vector se utilizar´a para generar el movimiento reactivo corrector, en caso de acercamiento; bien de la trayectoria o bien de la posici´on del robot. Para dise˜ nar el esquema de control, el primer paso es considerar la velocidad de acercamiento-alejamiento del obst´aculo (c) en un punto de la trayectoria (P ), dicha velocidad se define como la derivada temporal del gradiente: c=

d (∇u) dt

(4.13)

De esta velocidad se deduce una fuerza de colisi´on virtual de repulsi´on, dada por fV = K f · c (4.14)

donde Kf es la matriz de amortiguamiento, que relaciona el vector velocidad c con el vector fuerza de colisi´on virtual fV . Similarmente a como se hace en el caso de colisi´on f´ısica ([18] y [19]) interesa obtener la proyecci´on de esta fuerza sobre la superficie equipotencial que pasa por el punto P (representado por Seq en la figura 4.4) vr = kr (vT × (fV × vT )) ,

(4.15)

donde vT es la velocidad nominal de la trayectoria sobre una l´ınea de co∇u rriente (vT = vT |∇u| ) y kr es una constante real positiva. As´ı, el vector vr indica el camino de alejamiento de la trayectoria actual, situada sobre una l´ınea de corriente, a una nueva l´ınea de corriente que igualmente conducir´a al

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

71

Figura 4.5: Control correctivo por amortiguamiento. objetivo. El control s´olamente deber´a actuar en el caso de acercamiento del obst´aculo por debajo de una determinada distancia de seguridad dC . Es decir:   Constante > 0, para d · c < 0 y |d| < dC kr =  0, para d · c > 0 o |d| > dC

As´ı, la velocidad resultante de actuaci´on ser´a ve = vT − v r

(4.16)

Si esta vale cero (es decir el manipulador se halla en una posici´on deseada) y se le acerca un obst´aculo; este saldr´a de la posici´on actual con una velocidad proporcional a la velocidad de acercamiento. El sistema de control se representa en la figura 4.5, y tiene que cumplir que la velocidad de c´alculo de refrescamiento de datos debe ser superior a la velocidad de desplazamiento de los obst´aculos.

4.3.

Distribuci´ on espacial del cambio de la funci´ on potencial debido a un cambio local

En cualquiera de los casos: cambio de posici´on del punto destino, cambio de posici´on de un obst´aculo y aparici´on de un obst´aculo nuevo, la influencia que tendr´an sobre el cambio del valor del potencial ser´a progresivamente menos importante conforme mas alejados nos encontremos de la zona de cambio, y esta influencia ser´a tanto menor cuanto mayor sea el n´ umero de variables del C-espacio. Este hecho permite concentrar el c´alculo en las zonas y sus entornos pr´oximos donde se ha producido el cambio, reduci´endose el tiempo de c´omputo global.

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

72

Figura 4.6: Cambio de posici´on del punto destino. Cambio de posici´on del punto destino.- Suponiendo que se ha encontrado una trayectoria (partiendo de su funci´on potencial) con destino en el punto ri , y que se encuentra otra que enlace este punto con el situado en rj , interesa conocer la distribuci´on espacial del cambio de la funci´on potencial, debido u ´nicamente al cambio de posici´on del punto destino. Realizando la diferencia de modelos, se obtiene   ∇2 u = δ (r − ri ) en Ω  ∇2 u′ = δ (r − rj ) en Ω  (−) (=)   u = g en ∂Ω u = g en ∂Ω ∇2 v = δ (r − ri ) − δ (r − rj ) v = 0 en ∂Ω

 en Ω 

,



que es el modelo diferencia, con condiciones de contorno homog´eneas y donde v = u − u′ es la funci´on diferencia entre la nueva y la vieja funci´on potencial. Para evaluar, de forma cualitativa, la distribuci´on espacial del cambio de la funci´on potencial consideramos primero el entorno exterior a los dos puntos, el anterior y el actual, para ello aplicamos el teorema de la divergencia al modelo diferencial en una superficie que englobe a dichos puntos (figura 4.6) I Z Z
2 ~ ∇v·dS = ∇ v dV = (δ (r − ri ) − δ (r − rj )) dV = 0 (4.17) S

V

I

V

S

~= ∇v · dS

I

S

(∇v · n) dS =

I

S

∂v dS, ∂n

(4.18)

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

73

sustituyendo en la anterior I

S

∂v dS = 0 ∂n

(4.19)

y aplicando el teorema del valor medio para integrales   ∂v =0 ∂n M ED

(4.20)

que indica que el valor medio de la derivada direccional normal de la funci´on potencial diferencia, en torno de una superficie que engloba los dos puntos de cambio, es cero; considerando adem´as que las condiciones de contorno del modelo diferencial son homog´eneas, se tiene en cierto modo la condici´on cualitativa de que los efectos del cambio en entornos alejados es muy peque˜ na o nula. Tambi´en se puede ver que se cumple en cualquer punto del C-espacio (no coincidentes con los puntos ri y rj ) v = u − u′

⇒

∂v ∂u ∂u′ = − ∂n ∂n ∂n

∀n,

(4.21)

donde ∂/∂n es la derivada en la direcci´on n. Considerando el modelo ∇2 u = δ (r − ri ) con condiciones de contorno en el infinito, para una superficie de integraci´on Si centrada en el punto ri , se tiene   I 1 ∂u ∂u = (4.22) dS = 1 ⇒ ∂n M ED Si Si ∂n cuanto mayor es el n´ umero de variables del C-espacio y mayor es el radio de la hiperesfera de integraci´on m´as peque˜ no es el valor de la derivada direccional. El mismo razonamiento sirve para u′ y, utilizando la superficie de integraci´on S, se cumple I I I ∂v ∂u ∂u′ dS = dS − dS (4.23) S ∂n S ∂n S ∂n Con lo cual se concluye que el cambio de u a u′ se concentra principalmente (independientemente de las condiciones de contorno) en los entornos pr´oximos de las zonas donde ha desaparecido δ(ri ) y donde ha aparecido δ(rj ). Cambio de posici´on de un obst´aculo.- En este caso, a diferencia del apartado 4.1, no cuenta la variable tiempo; mientras dura el proceso de c´alculo de una trayectoria, las posiciones tanto del punto destino como de los obst´aculos se consideran invariables, constituyendo un escenario de c´alculo (i). La posibilidad de que alg´ un obst´aculo cambie (normalmente de forma ligera) su posici´on en el C-espacio da lugar a un escenario de c´alculo nuevo (i + 1).

74

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

Figura 4.7: Cambio de posici´on de un obst´aculo. En ambas posibilidades de modelado de los obst´aculos, como contorno y como funci´on, se obtiene el mismo resultado: la influencia del cambio se aten´ ua con la distancia, tanto m´as cuanto mayor es el n´ umero de grados de libertad. En la figura 4.7, el obst´aculo O se desplaza ligeramente de la posici´on Oi a la posici´on Oi+1 de forma que una zona ser´a coincidente y otra desaparece de la fase i para aparecer en la fase i + 1, u ´nicamente esta zona de cambio ser´a la que afecte a la variaci´on de la funci´on potencial en el nuevo escenario de c´alculo. - Obst´ aculo modelado como contorno interior.- Para un obst´aculo modelado como contorno interior a un dominio Ω  ∇2 u = δ (r − rd ) en Ω  ,  u = f en ∂Ω

donde el obst´aculo O est´a definido por el contorno ∂ΩO . La soluci´on ser´a u(r) = g(r, rd ) + IO (4.24)

La aportaci´on (como suma) a la soluci´on dada por el obst´aculo viene dada por la integral I ∂G dS, (4.25) f (rc ) IO = ∂n ∂ΩO donde G es la funci´on de Green y rc representa la posici´on de un punto cualquiera del contorno del obst´aculo. Si f (rc ) = k, ∀rc ∈ ∂ΩO , la integral valdr´a I ∂G IO = k dS = k · 1 = k, ∀r ∈ / O, (4.26) ∂ΩO ∂n

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

75

esto es as´ı independientemente de la posici´on del obst´aculo. En cuanto al otro t´ermino de la soluci´on g(r − rd ), es coincidente con la funci´on de Green situada en rd y, atendiendo a su m´etodo de obtenci´on mediante las cargas imagen, existir´a una carga imagen neta negativa (−δ) en el interior del obst´aculo, de forma que aplicando el teorema de la divergencia a una superficie que contenga dicho obst´aculo se cumplir´a I Z ∇u · dS = (−δ) dV = −1, (4.27) S

I

S

V

∇u · dS =

I

S

(∇u · n) dS =

I

S

∂u dS = −1 ∂n

aplicando el teorema del valor medio   1 ∂u =− ∂n M ED S

(4.28)

(4.29)

que indica que el valor medio de la derivada direccional normal de la funci´on potencial, en torno de una superficie que engloba al obst´aculo, es inversamente proporcional a dicha superficie. En conclusi´on: el efecto que tiene en la funci´on potencial un cambio de posici´on de un obst´aculo (modelado como contorno interior) es apreciable en el entorno pr´oximo de (figura 4.6) Oǫ = (Oi ∪ Oi+1 ) − (Oi ∩ Oi+1 ) ,

(4.30)

y es muy peque˜ no en el entorno situado a partir de cierta distancia. - Obst´ aculo modelado como funci´on.- Para el caso del movimiento de un obst´aculo modelado como funci´on, se obtiene el modelo diferencia   ∇2 u′ = ρi+1 (r) en Ω  ∇2 u = ρi (r) en Ω  (=) (−)   u = g en ∂Ω u = g en ∂Ω ∇2 v = ρǫ (r)

 en Ω 

v = 0 en ∂Ω donde ρǫ =

  > 0, 

= 0,

,



∀r ∈ Oǫ ∀r ∈ / Oǫ

Aplicando el teorema de la divergencia al modelo diferencia, con una

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

76

superficie que englobe al obst´aculo funci´on en las dos posiciones (figura 4.6) I Z S

∇v · dS =

ρǫ (r) dV

(4.31)

V

Suponiendo que la funci´on ρ que modela el obst´aculo tiene un valor constante (ρ = k)   I ∂v ∂v kVǫ (r) dS = kVǫ (r) ⇒ , (4.32) = ∂n M ED S S ∂n

donde Vǫ es el volumen diferencial del movimiento (Oǫ ) y S es la superficie de integraci´on, cuyo radio es considerado como variable para indicar la decreciente influencia del cambio conforme nos alejamos del mismo. Aparici´on de un nuevo obst´aculo.- La exploraci´on continuada del Cespacio puede dar lugar a la aparici´on de nuevos obst´aculos y, por tanto, al cambio del entorno de c´alculo. La aparci´on de un obst´aculo en el entorno tiene un tratamiento y unas consecuencias similares al movimiento de un obst´aculo; la influencia que tendr´a en la variaci´on de la funci´on potencial ser´a importante en el entorno del mismo y ser´a cada vez m´as d´ebil conforme nos alejamos del mismo. Para el caso de un obst´aculo modelado como contorno interior el efecto sobre la funci´on potencial ser´a, en cuanto a la integral de contorno I ∂G f (rc ) IO = dS = (f (rc ))M ED , (4.33) ∂n ∂ΩO y en cuanto a su aportaci´on a la funci´on de Green se cumple igualmente   1 ∂u (4.34) =− , ∂n M ED S donde S es cualquier superficie que rodee el obst´aculo, por ejemplo una funci´on esf´erica de radio R. Obst´ aculo modelado como funci´on.- Si un determinado entorno, modelado por ∇2 u = ρ (r) y condiciones de contorno en el infinito, experimenta una variaci´on con la aparici´on de un obst´aculo modelado por la funci´on ρO , el valor de la funci´on potencial cambiar´a ∇2 u′ = ρ (r) + ρO (r)

(4.35)

Para ver la influencia de este cambio a una cierta distancia del nuevo obst´aculo se restan las dos ecuaciones anteriores ∇2 (u′ − u) = ρO (r) ,

(4.36)

y aplicando el teorema de la divergencia en una superficie esf´erica, de radio R, que contenga a ρO I Z Z
2 ∇v · dS = ∇ v dV = ρO (r) dV = ρO VO , (4.37) S

V

V

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

77

para v = u′ − u. Por la simetr´ıa del modelo planteado, el gradiente ser´a de m´odulo constante y normal a la superficie esf´erica de integraci´on en cada punto de la misma. (|∇v|)M ED =

ρO V O S

(4.38)

El m´odulo del gradiente, debido al nuevo obst´aculo aparecido, es inversamente proporcional a Rn−1 para n dimensiones. Por lo tanto, la influencia del nuevo obst´aculo en la funci´on potencial ser´a muy peque˜ na a partir de una cierta distancia, con lo cual s´olamente ser´a necesario refrescar el c´alculo en una peque˜ na zona en torno del mismo, reduci´endose el tiempo de c´omputo.

4.4.

Aportaci´ on

En este cap´ıtulo se ha particularizado el estudio delas Funciones Arm´onicas para la planificaci´on de movimientos en entornos din´amicos, llegando a las siguientes conclusiones: - Cuando los entornos son din´amicos y los obst´aculos pueden moverse a una determinada velocidad, es posible obtener una funci´on potencial superarm´onica (sin m´ınimos locales) modelando el obst´aculo como una funci´on generadora de calor que se desplaza a una determinada velocidad; dimensionando adecuadamente la constante de difusi´on, es posible que la derivada temporal venga acotada, que tengamos una funci´on potencial din´amica pero siempre superarm´onica, de forma que no presente m´ınimos locales. Adem´as, utilizando esta funci´on potencial din´amica, se propone un esquema de control reactivo virtual que sirva para corregir las trayectorias (y configuraciones est´aticas de reposo) evitando las colisiones f´ısicas. - Cuando los entornos son cambiantes; bien por el cambio de posici´on del punto destino, bien por el desplazamiento de un obst´aculo o por la aparici´on de uno nuevo, el cambio correspondiente de la funci´on potencial se localiza principalmente en el entorno pr´oximo de aquellas zonas que han cambiado y a partir de una cierta distancia la percepci´on de este cambio es m´as d´ebil. El tiempo de computaci´on global se ver´a reducido si se dedica un mayor n´ umero de iteraciones de c´omputo a refrescar el valor del potencial en el entorno pr´oximo donde se ha producido el cambio.

78

´ EN ENTORNOS DINAMICOS ´ CAP. 4. PLANIFICACION

Cap´ıtulo 5 C´ alculo num´ erico de las funciones arm´ onicas, superarm´ onicas y subarm´ onicas Para obtener la soluci´on num´erica de la ecuaci´on de Poison se deduce una ecuaci´on en diferencias definida sobre una discretizaci´on del espacio. Para la resoluci´on de esta ecuaci´on en diferencias existen dos m´etodos: I) M´etodos directos, que utilizan una notaci´on matricial, con resoluciones algebraicas como son los m´etodos de Gauss y sus derivados. Est´an limitados en cuanto al n´ umero de puntos de c´alculo utilizados. II) Los m´etodos de relajaci´on de Jacobi, de Gauss-Seidel y de SOR, se usan cuando el n´ umero de puntos de c´alculo es muy elevado ya que permiten obtener una soluci´on v´alida con un n´ umero de iteraciones limitado. Se corresponden con la resoluci´on cl´asica de ecuaciones en diferencias lo que ofrece una visi´on f´ısica del m´etodo, permitiendo extraer expresiones anal´ıticas sobre las que se deducen propiedades utilizables en la reducci´on del tiempo de c´alculo. La multirresoluci´on es un m´etodo donde la ecuaci´on en diferencias, utilizada para el c´alculo, se define sobre rejillas uniformes con distintas resoluciones; el c´alculo coordinado y correctivo, en las distintas rejillas, permite reducir el tiempo de c´omputo bas´andose en las propiedades de los m´etodos de relajaci´on. En este cap´ıtulo se propone el c´alculo de la ecuaci´on en diferencias sobre una discretizaci´on jer´arquica en multirresoluci´on del espacio, lo que permite una reducci´on importante del n´ umero de puntos de c´alculo y el aprovechamiento de las propiedades de la multirresoluci´on. En primer lugar, en el apartado 5.1, se obtiene la expresi´on de la ecuaci´on en diferencias, sobre una discretizaci´on uniforme; que servir´a de base para la definici´on de los distintos m´etodos de relajaci´on. 79

´ ´ CAP. 5. CALCULO NUMERICO

80

Despu´es, en el apartado 5.2, se introduce como nuevo m´etodo de c´alculo, la discretizaci´on y multirresoluci´on jer´arquica adaptativa.

5.1.

Discretizaci´ on uniforme

En una discretizaci´on uniforme el espacio queda particionado por una rejilla constituida por celdillas regulares c de lado h, en las que la funci´on potencial toma un determinado valor.

5.1.1.

Ecuaci´ on en diferencias

A una determinada rejilla regular de resoluci´on h (Ωh ), le corresponde una plantilla de c´alculo mostrada en la figura 5.1. Utilizando Taylor para dos variables, el valor del potencial en un punto (q1 , q2 ) respecto del potencial en otro punto pr´oximo (q10 , q20 ), es: 1 ∂u 1 ∂u (q1 − q10 ) + (q2 − q20 )+ u(q1 , q2 ) = u(q10 , q20 ) + 1! ∂q1 (q10 ,q20 ) 1! ∂q2 (q10 ,q20 ) 1 ∂ 2 u 1 ∂ 2 u 2 (q1 − q10 ) + (q2 − q20 )2 + 2! ∂q12 (q10 ,q20 ) 2! ∂q22 (q10 ,q20 ) 2 ∂ 2 u (q1 − q10 )(q2 − q20 ) + . . . (5.1) 2! ∂q1 ∂q2 (q10 ,q20 )

Aplicado a la disposici´on en cruz de la figura 5.1:

 2 3 ∂u u(q10 + h, q20 ) = u(q10 , q20 ) + 1!1 ∂q h + 2!1 ∂∂qu2 h2 + 3!1 ∂∂qu3 h3 + . . .   1 1 1       1 ∂2u 2 1 ∂u 1 ∂3u 3    u(q10 − h, q20 ) = u(q10 , q20 ) − 1! ∂q1 h + 2! ∂q2 h − 3! ∂q3 h + . . . 1

1

 2 3 ∂u   h + 2!1 ∂∂qu2 h2 + 3!1 ∂∂qu3 h3 + . . . u(q10 , q20 + h) = u(q10 , q20 ) + 1!1 ∂q  2  2 2       u(q , q − h) = u(q , q ) − 1 ∂u h + 1 ∂ 2 u h2 − 1 ∂ 3 u h3 + . . . 10 20 10 20 1! ∂q2 2! ∂q 2 3! ∂q 3 2

2

donde las derivadas est´an evaluadas en el punto (q10 , q20 ). Sumando las cuatro ecuaciones, tenemos la siguiente ecuaci´on en diferencias: u(q10 + h, q20 ) + u(q10 − h, q20 ) + u(q10 , q20 + h) + u(q10 , q20 − h) =   2 ∂ 2u ∂ 2u 4u(q10 , q20 ) + + 2 h2 (5.2) 2! ∂q12 ∂q2 ya que los t´erminos correspondientes a las derivadas de orden impar se anulan, y se consideran nulas la suma de derivadas de orden cuatro y superiores.

´ ´ CAP. 5. CALCULO NUMERICO

81

Figura 5.1: Discretizaci´on uniforme: Disposici´on 4 vecinos en cruz. Generalizando para n dimensiones y considerando la discretizaci´on del espacio con N divisiones por cada variable: qj = (j1 , j2 , . . . , jl , . . . , jn )h | jl ∈ {0, 1, . . . , N − 1} 2n X

uqji = 2nuqj +

i=1


2h2 ∇2 u q , j 2!

(5.3) (5.4)

representando con el sub´ındice i cada una de las celdillas vecinas de acuerdo con las disposiciones en cruz, obteni´endose 2n

uq j


1 X h2 = ∇2 u q uqji − j 2n i=1 2!n

∀qj ∈ Ωh

(5.5)

Siendo esta ecuaci´on en diferencias utilizable computacionalmente en el c´alculo de la funci´on en un espacio, el cual dispondr´a de unas determinadas condiciones de contorno. Para las funciones arm´onicas se cumple la ecuaci´on de Laplace (∇2 u = 0), v´alida en todo punto del dominio, obteni´endose 2n

uq j

1 X uq , = 2n i=1 ji

(5.6)

en donde se cumple la propiedad del valor medio, por consiguiente, en ning´ un punto del dominio existir´a un m´ınimo, ni un m´aximo. La ecuaci´on de Poison (∇2 u = ρ(r)) modela las funciones superarm´onicas (ρ(r) < 0) y la subarm´onicas (ρ(r) > 0). As´ı, se obtiene 2n

uq j donde: ρ0 = ρ(qj )

1 X h2 = uqji − ρ0 , 2n i=1 2!n

(5.7)

´ ´ CAP. 5. CALCULO NUMERICO

82

Figura 5.2: C´alculo del potencial de una celda vecina a otra obst´aculo, en una discretizaci´on uniforme. (b) Condiciones de contorno de Dirichlet. (c) Condiciones de contorno de Neuman.

5.1.2.

Condiciones de contorno de Dirichlet y de Neuman

Para condiciones de contorno de Dirichlet, dadas por u(r) = uC ,

∀r ∈ ∂Ω,

(5.8)

aplicadas en un espacio discretizado uniformemente (figura 5.2a), las celdas correspondientes a la discretizaci´on del contorno afectan al valor del potencial de las celdas vecinas del espacio libre, de forma que un n´ umero determinado de celdillas dentro del sumatorio pertenecer´an al contorno de C-obst´aculo, cumpli´endose que ! NL NC X 1 X h2 ρ0 , (5.9) uq j = uqji + uC − 2n i=1 2!n j=1 donde uqj (u0 en la figura 5.2b) es el potencial de una celdilla vecina del contorno, uqji es el potencial de una celdilla vecina libre y uC es el potencial de la celdilla vecina contorno. NL es el n´ umero de celdillas vecinas libres y NC el de celdillas contorno NL + NC = 2n

(5.10)

Para condiciones de contorno de Neuman, dadas por ∂u (r) = Nu , ∂n

∀r ∈ ∂Ω,

(5.11)

aplicadas en un espacio discretizado uniformemente, el vector normal al contorno (n) de un C-obst´aculo discretizado, se refiere a las celdillas que contienen el contorno, de manera que en un punto queda adjudicado a la celdilla

´ ´ CAP. 5. CALCULO NUMERICO

83

que lo contiene y apunta a una u ´nica celdilla contigua del espacio libre (figura 5.2c), a la que influye en su c´alculo de la manera siguiente: uobs − uqj = Nu , ∆q

(5.12)

uobs = uqj + Nu ∆q

(5.13)

de donde y uqj

1 = NL

NL X i=1

uqji +

NC X i=1

uobs

!

1 = NL

NL X i=1

uqji +

NC X i=1

Nu ∆q

!

(5.14)

Donde ∆q = h. Lo habitual es que Nu valga cero, en cuyo caso el c´alculo se simplifica: las celdillas de contorno quedan excluidas del c´alculo del valor medio, en sus celdillas vecinas del espacio libre. En contrapartida las trayectorias se aproximar´an a los obst´aculos.

5.1.3.

M´ etodos iterativos

Son los m´etodos num´ericos utilizados para encontrar la soluci´on de la ecuaci´on de Poison y se corresponden, b´asicamente, con la discretizaci´on de la ecuaci´on del calor: ∂u = ∇2 u − ρ(r) (5.15) ∂t Aplicando diferencias progresivas en el tiempo, para la obtenci´on de la derivada temporal k uk+1 ∂u qj − uqj ≈ (5.16) ∂t ∆t Utilizando las ecuaciones (5.5), (5.15) y (5.16), se obtiene el modelo gen´erico de los m´etodos iterativos # " 2n X k ukqji − 2nukqj − (∆t) ρ0 , ∀qj ∈ Ωh uk+1 (5.17) q j = uq j + s i=1

donde s=

∆t h2

(5.18)

y n es el n´ umero de variables. En r´egimen permanente no hay variaci´on temporal de la funci´on k uk+1 q j = uq j = uq j

(5.19)

y la soluci´on de la ecuaci´on del calor se corresponde con la soluci´on de la ecuaci´on de Poison expresada en (5.7)

´ ´ CAP. 5. CALCULO NUMERICO

84

Figura 5.3: Eigenvalores e iteraciones en funci´on de n´ umero de onda para el m´etodo de Jacobi. 5.1.3.1.

M´ etodo iterativo de Jacobi

Haciendo (s = 1/2n) y, teniendo en cuenta que el tiempo es irrelevante en este caso en cuanto a que no afecta a la soluci´on, (∆t = 1) en la ecuaci´on (5.17) se obtiene el m´etodo recursivo de Jacobi para la ecuaci´on de Poison (para la ecuaci´on de Laplace ρ0 vale cero) 2n

ukqj =

1 X k−1 − ρ0 , u 2n i=1 qji

∀qj ∈ Ωh

(5.20)

obteni´endose valores cada vez m´as aproximados a la soluci´on buscada. Restando la ecuaci´on aproximada de la exacta: P2n  1 uqj = 2n  i=1 uqji − ρ0 , P2n k−1  1 k uqj = 2n i=1 uqji − ρ0

se obtiene la ecuaci´on de diferencias del error: 2n

ekqj

1 X k−1 e , = 2n i=1 qji

(5.21)

donde (ekqj = uqj − ukqj ) tiende a cero conforme ukqj tiende a uqj y 2n

uq j

1 X k−1 − ρ0 + ekqj u = 2n i=1 qji

(5.22)

La soluci´on de la ecuaci´on (5.22) se obtiene en el ap´endice (A.5.1) y vale ekqj

n N −1 Y X

p π  l λ sin = jl N p =1 l=1 l

k

(5.23)

´ ´ CAP. 5. CALCULO NUMERICO

85

donde N es el n´ umero de particiones de la discretizaci´on, pl es el modo o n´ umero de onda (pl = 1, 2, 3, . . . , N − 1) y λ son los eigenvalores (figura 5.3) n

p π  1X l λ= cos n l=1 N

(5.24)

La soluci´on del error se muestra como una serie de Fourier (para una funci´on de n variables) cuyas componentes se van aproximando iterativamente a cero. Puesto que (−1 < λ < 1), se cumple que (l´ımk→∞ e = 0) y el m´etodo converge (ap´endice A.5.1 ecuaci´on A.110), a una velocidad determinada por: log 2 , (5.25) π2 donde K es el n´ umero de iteraciones necesarias para reducir el error a la mitad. Se aprecia que la velocidad de convergencia es muy baja, independiente del n´ umero de variables. De la expresi´on anterior tambi´en se deduce que cuanto menor sea N (rejilla mas gruesa) mayor ser´a la velocidad de convergencia. Esta cuesti´on es muy importante a la hora de dise˜ nar algoritmos en multirresoluci´on. Adem´as, si los arm´onicos de la ecuaci´on del error (5.24) tienen un n´ umero de onda (pl ) pr´oximo a N/2, la atenuaci´on de los mismos es m´as r´apida (figura 5.3) K = 2N 2

5.1.3.2.

M´ etodo iterativo de Gaus-Seidel

Este m´etodo modifica el de Jacobi de forma que utiliza los valores de la funci´on en aquellas celdillas vecinas, que justamente se acaban de calcular en la iteraci´on presente, en lugar de utilizar los valores correspondientes a la iteraci´on anterior. Con este fin el orden de c´alculo, de celdillas, se realiza en diagonal: se van incrementando correlativamente cada una de las variables. As´ı, este m´etodo utiliza, para la ecuaci´on de Poisson, la expresi´on de c´alculo siguiente: n n 1 X k−1 1 X k k + − ρ0 (5.26) u u uqj = 2n i=1 qji− 2n i=1 qji+ donde ukqji− representa los valores de las celdillas vecinas, que se han calculado en la iteraci´on actual k y uk−1 qji+ representa los valores de las celdillas vecinas, que se han calculado en la iteraci´on anterior (k − 1). Entonces, si se resta (como antes) la expresi´on recursiva de la exacta  P2n 1 uqj = 2n  i=1 uqji − ρ0 , Pn k Pn k−1  1 1 k uqj = 2n i=1 uqji− + 2n i=1 uqji+ − ρ0 se obtiene

 k  eqj = 

uq j =

1 2n 1 2n

Pn

k i=1 eqji−

Pn

i=1

+

ukqji− +

1 2n 1 2n

Pn

k−1 i=1 eqji+

Pn

i=1

k uk−1 qji+ − ρ0 + eqj

´ ´ CAP. 5. CALCULO NUMERICO

86

Figura 5.4: Eigenvalores en funci´on de n´ umero de onda para el m´etodo GaussSeidel y N = 50.

Tal como se muestra en el ap´endice (A.5.2), la soluci´on es como antes (ecuaci´on 5.24), pero los eigenvalores, en este caso, vienen dados por |λ| =

s

1 , pl = 1, 2, 3, . . . , N − 1 , 5 − 4 cos( pNl π )

(5.27)

sus valores, en funci´on del n´ umero de onda, vienen representados en la figura 5.4. Tambi´en en este caso (−1 < λ < 1) y (l´ımk→∞ e = 0); en el ap´endice (A.5.2), ecuaci´on (A.129), se demuestra que el n´ umero de iteraciones necesario para reducir el error a la mitad vale K = N2

log 2 , π2

(5.28)

lo que supone una reducci´on respecto al m´etodo de Jacobi de 2. En este caso, como se aprecia en la figura 5.4, cuanto mayor es el n´ umero de onda, menor es el valor de λ y m´as r´apidamente se aten´ ua el arm´onico correspondiente; actuando como un filtro paso bajo. As´ı, si se establecen unas condiciones iniciales aleatorias, con componentes en alta frecuancia, m´as r´apidamente se alcanzar´a el valor deseado. 5.1.3.3.

M´ etodo iterativo S.O.R. (Sucesive Over Relaxation)

Este m´etodo incorpora, respecto del anterior, un par´ametro de relajaci´on (ω) de la manera siguente: " n # n X X ω ukqj = (1 − ω)uk−1 uk + uk−1 − ρ0 (5.29) qj + 2n i=1 qji− i=1 qji+

´ ´ CAP. 5. CALCULO NUMERICO

87

Figura 5.5: Eigenvalores en funci´on de n´ umero de onda para el m´etodo S.O.R. y N = 50.

para la ecuaci´on de Laplace ρ0 vale cero. Nuevamente, restando la de relajaci´on de la expresi´on exacta hP i Pn n ω uqj = (1 − ω)uqj + 2n i=1 uqji− + i=1 uqji+ − ρ0 ukqj

= (1 −

ω)uk−1 qj

+

ω 2n

se obtiene uqj = ukqj + ekqj

hP

n i=1

ukqji−

+

Pn

k−1 i=1 uqji+

i

− ρ0

expresi´on    

,

  

" n # n X X ω = (1 − ω)uk−1 uk + uk−1 − ρ0 + ekqj (5.30) qj + 2n i=1 qji− i=1 qji+

y la expresi´on del error n

ekqj

= (1 −

ω)ek−1 qj

n

ω X k−1 ω X k eqji− + e + 2n i=1 2n i=1 qji+

(5.31)

cuyos eigenvalores (de la soluci´on se obtienen en el ap´endice (A.5.3)) valen |λ| =

s

(1 − ω)2 + ( ω2 )2 + (1 − ω)ω cos( pNl π ) 1 + ( ω2 )2 − ω cos( pNl π )

(5.32)

En la figura 5.5 se muestra los eigenvalores en funci´on del n´ umero de onda, para distintos valores del par´ametro ω. Se observa, tambi´en, como los arm´onicos, componentes del error, de alta frecuencia se atenuan m´as r´apidamente, existiendo un valor ´optimo de ω que minimiza |λ|2 : r    √ pl π 1 (5.33) 2 − 2 2 1 − cos ω= 2cos pNl π − 1 N

88

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.6: Aumento de la resoluci´on de la rejilla e interpolaci´on de puntos. El c´alculo de los nuevos nodos, en rojo, se realiza en dos fases: I Vecinos en X, figura (b), II Vecinos en cruz, figura (c). Para (N = 50) y (pl = 1), su valor es 1,88. Y, el n´ umero de iteraciones que reduce a la mitad el error (A.5.3), ecuaci´on (A.143) viene dada por   2 log 2 K=N (5.34) π que es un n´ umero mucho m´as reducido que los anteriores (ecuaciones 5.26 y 5.29 ) ya que N no viene elevado al cuadrado. Por este motivo, el m´etodo iterativo S.O.R. es el utilizado en el planificador de trayectorias.

5.1.4.

C´ alculo en multirresoluci´ on

Puesto que λ es la amplitud de los arm´onicos de Fourier que componen el error y es, al mismo tiempo, la base de la expresi´on recursiva temporal, los m´etodos iterativos act´ uan sobre dichos arm´onicos de manera desigual. Excepto el de Jacobi, amortiguan m´as r´apidamente los de alta que los de baja frecuencia. Por otro lado las bajas frecuencias del error en una determinada rejilla de discretizaci´on fina, son frecuencias altas en un discretizaci´on gruesa. De esta manera se pueden realizar inicialmente iteraciones sobre una discretizaci´on gruesa, amortiguando los arm´onicos m´as r´apidamente, y pasar despu´es los valores obtenidos a una discretizaci´on fina sobre la que se completen los c´alculos iterativos. Adem´as el n´ umero de puntos de c´alculo se reduce en un factor de 1/2kn , donde n es el n´ umero de variables y k el n´ umero de veces que una divisi´on gruesa es mayor que la fina. Al pasar de una discretizaci´on gruesa a otra m´as fina aparecen nuevos puntos en los cuales el valor de la funci´on es desconocido. Entonces, se realiza una estimaci´on inicial en dichos puntos realizando una interpolaci´on (figuras 5.6a, b y c) lineal y a partir de aqu´ı se sigue calculando sobre la nueva rejilla. De la misma forma que el valor exacto del potencial es desconocido, tambi´en lo es el error; entonces para relacionar el error con la funci´on, en cada

´ ´ CAP. 5. CALCULO NUMERICO

89

Figura 5.7: Correcci´on con rejilla gruesa. iteraci´on, se define el residuo (en una celdilla i) como: rik = uki − uk−1 i

| k = 1, 2, 3, . . .

(5.35)

Para obtener una relaci´on iterativa del error, relacionada con el residuo, se tiene en cuenta que: eki = ui − uki (5.36) y ) = −eki + ek−1 rik = (ui − eki ) − (ui − ek−1 i i

(5.37)

de donde, se obtiene: eki = ek−1 − rik i

(5.38)

para cada celdilla i. La idea b´asica del c´alculo en multirresoluci´on es la siguiente (figura 5.7): 1. Realizar unas pocas iteraciones sobre una rejilla de resoluci´on h (Ωh ) para obtener una aproximaci´on de la funci´on u y del error e. Por ejemplo, para el m´etodo de Jacobi se utilizar´ıan las ecuaciones (5.20),(5.21) y (5.22). Donde los valores iniciales u0i se establecen de forma aleatoria, ya que su espectro de arm´onicos de Fourier es uniforme y, por tanto, sus componentes de alta frecuencia ser´an atenuados r´apidamente con las primeras iteraciones. 2. Disminuir la resoluci´on a 2h y transferir (mapear) el error: e(2h) ← e(h) 3. Realizar unas pocas iteraciones sobre Ω2h , utilizando la ecuaci´on (5.21), para obtener una aproximaci´on del error e. 4. Aumentar la resoluci´on, nuevamente, a h e interpoler los valores del error: Ω2h → Ωh 5. Corregir la aproximaci´on: uk ← uk + ek De una manera recursiva se puede ir relajando y mapeando de rejillas m´as finas a rejillas m´as gruesas, para seguidamente volver interpolando y corrigiendo siguiendo un ciclo en V (figura 5.8). Se consigue una soluci´on con

90

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.8: Multirresoluci´on: (a) Ciclo V. (b) Ciclo FMG. menos iteraciones si los valores iniciales son buenos, para ello se comienza relajando en la resoluci´on m´as gruesa y sucesivamente se va interpolando hacia rejillas m´as finas donde se vuelve a relajar, se mapea y vuelve a relajar en la rejilla gruesa. Sucesivamente se va aumentando la amplitud del ciclo hasta alcanzar la soluci´on en la rejilla fina. Este m´etodo recursivo se denomina Full Multigrid (FMG) [15].

5.2.

Discretizaci´ on jer´ arquica adaptativa

Puesto que el objetivo es obtener una funci´on potencial auxiliar que sirva para obtener una trayectoria por seguimiento del gradiente, no se precisa resoluci´on excepto en los contornos de los obst´aculos. Por este motivo resulta ´optima una discretizaci´on jer´arquica adaptativa; la cual, sucesivamente, subdivide las celdas que contienen contorno de obst´aculo hasta alcanzar la m´axima resoluci´on prefijada (figura 5.9). Con esta discretizaci´on se obtiene una reducci´on m´axima de los puntos de c´alculo ya que en los espacios vac´ıos se adaptar´an celdas de m´aximo tama˜ no, representando un solo punto de c´alculo. Las ventajas obtenidas con esta discretizaci´on son: • Optimiza la discretizaci´on fina, concentr´andola en el contorno de los obst´aculos. • Con una resoluci´on m´axima en los contornos, puntos inicial y final, el n´ umero de puntos de c´alculo se reduce en los espacios vac´ıos y s´olidos. • El n´ umero de iteraciones necesarias en los m´etodos de relajaci´on es menor cuanto menor es el n´ umero de celdas de c´alculo. Para obtener esta discretizaci´on recursiva y adaptativa se precisa una funci´on detectora de colisiones que reporte la distancia de colisi´on; con ella se

´ ´ CAP. 5. CALCULO NUMERICO

91

Figura 5.9: Niveles de discretizaci´on jer´arquica adaptativa al contorno.

Figura 5.10: C debe identificar univocamente una celda, en tama˜ no y posici´on. determina si una celda corresponde a una celda libre, obst´aculo o mixta.

5.2.1.

Esquema num´ erico de identificaci´ on de celdas

Para poder calcular el valor de la funci´on en cada celda de la discretizaci´on, es preciso ubicar su posici´on en el espacio. Mientras que en una discretizaci´on regular esta queda determinada por sus coordenadas, en una discretizaci´on jer´arquica es necesario establecer un convenio num´erico que identifique un´ıvocamente una celda, tanto su tama˜ no como su situaci´on en el espacio de configuraciones. As´ı, en la figura 5.10, a una determinada celda c le corresponder´an unas coordenadas (q1 , q2 ) y un determinado tama˜ no. Con la finalidad de poder representar cada celda, se establecen distintos niveles de partici´on (m ∈ N) tal que (0 ≤ m ≤ M ), donde M ser´a el m´aximo nivel de partici´on.

92

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.11: Rejillas para diferentes niveles de partici´on, en un espacio de configuraciones 2D. Las celdas marcadas se refieren al ejemplo expresado en la ecuaci´on (5.48) Y, para un espacio de n dimensiones, a cada nivel de partici´on m le corresponder´a n una rejilla regular denominada Gm (figura 5.11), que contendr´a 2nm celdas del tipo cnm , con sus coordenadas cartesianas: qm = (qm1 , qm2 , . . . , qmi , . . . , qmn )

(5.39)

De esta forma, cualquier celdilla de una discretizaci´on jer´arquica, vendr´a univocamente identificada por su nivel de partici´on m (tama˜ no) y sus coorden nadas dentro de la rejilla Gm , a la que pertenece: c = cnm ≡ {m, (qm1 , qm2 , . . . , qmi , . . . , qmn )}

(5.40)

Y, el conjunto de rejillas necesarias para la discretizaci´on jer´arquica, en un espacio de n dimensiones con un nivel m´aximo de partici´on M , ser´a: n G n = (G0n , G1n , . . . , GM ),

(5.41)

Por ejemplo, en la figura 5.10 la celda c pertenece a la rejilla G32 y sus coordenadas son (3, 4), es decir: c = c23 ≡ {m, (qm1 , qm2 )} = {3, (3, 4)}

(5.42)

En una discretizaci´on jer´arquica adaptativa, si una celda cnm debe subdividirse, generar´a 2n celdas del tipo cnm+1 (incluidas en ella), cuyas coordenadas valdr´an: q(m+1)i = 2 · qmi + l(m+1)i | l(m+1)i ∈ {0, 1}, (5.43)

∀i ∈ {1, 2, . . . , n} y m ∈ {0, 1, 2, . . . , M − 1} Inversamente, cualquier celda cnm estar´a incluida en otra del tipo cn(m−1) de coordenadas: q  mi q(m−1)i = int ∀i ∈ {1, 2, . . . , n} ∀m ∈ {0, 1, 2, . . . , M − 1} (5.44) 2

´ ´ CAP. 5. CALCULO NUMERICO

93

De esta forma, la simplicidad de codificaci´on de celdas presentada, junto con la utilizaci´on de la Multirresoluci´on Jer´arquica Adaptativa, permitir´a una buena eficiencia de computaci´on de la ecuaci´on de Poison. Las coordenadas de las celdas del tipo cn(m+h) , contenidas en una determinada celda cnm , se obtienen utilizando recursivamente la ecuaci´on (5.47): q(m+h)i = 2h qmi +2h−1 l(m+1)i +2h−2 l(m+2)i +. . .+21 l(m+h−1)i +20 l(m+h)i (5.45) y, expresado matricialmente:      l(m+1)1 l(m+2)1 qm1 q(m+h)1  qm2   l(m+1)2 l(m+2)2  q(m+h)2        = 2h  .. +  .. .. ..   .    . . . l(m+1)n l(m+2)n qmn q(m+h)n

. . . l(m+h)1 . . . l(m+h)2 .. ... . . . . l(m+h)n

    

2(h−1) 2(h−2) .. .

    

20 (5.46)

donde

Anm+h



  = 

l(m+1)1 l(m+2)1 l(m+1)2 l(m+2)2 .. .. . . l(m+1)n l(m+2)n

. . . l(m+h)1 . . . l(m+h)2 .. ... . . . . l(m+h)n

    

(5.47)

es la matriz binaria de coordenadas, cuyas columnas representan las coordenadas de ubicaci´on en cada nivel con respecto al anterior. Por ejemplo, en la figura 5.11, la celda del tipo c23 de coordenadas (6, 1), est´a contenida en otra del tipo c21 de coordenadas (1, 0) y relacionadas entre ellas por la expresi´on:       1    q(3)1 6 1 1 0 2 2 = (5.48) =2 + 0 q(3)2 1 2 0 0 1 donde la matriz binaria de coordenadas   1 0 2 A2+1 = 0 1

(5.49)

tiene en su segunda columna las coordenadas (0, 1) correspondientes a la posici´on de la celda (6, 1) del nivel 3, respecto de la celda (3, 0) del nivel 2 que la contiene y en la primera columna tiene las coordenadas (1, 0) correspondientes a la posici´on de la celda (3, 0) del nivel 2, respecto de la celda (1, 0) del nivel 1 que la contiene. Rosell e I˜ niguez, en la publicaci´on [1] del apartado 7.1, proponen un m´etodo de etiquetado y c´alculo de las funciones arm´onicas, utilizadas en la planificaci´on de movimientos en rob´otica.

5.3.

Multirresoluci´ on Jer´ arquica Adaptativa

En esta tesis se propone utilizar la t´ecnica de multirresoluci´on aplicada a la discretizaci´on jer´arquica. Para ello es necesario definir un nuevo m´etodo de multirresoluci´on denominado Multirresoluci´on Jer´arquica Adaptativa, que se basar´a en los principios siguientes:

94

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.12: Multirresoluci´on Jer´arquica Adaptativa. Las figuras de la columna (a) representan las ditintas fases (de arriba a abajo) de dicretizaci´on jer´arquica, y las de las columna (b) son sus correspondientes para el c´alculo en discretizaci´on uniforme.

´ ´ CAP. 5. CALCULO NUMERICO

95

Tal como se muestra en la figura 5.12, las rejillas de la multirresoluci´on uniforme se adaptan a la discretizaci´on jer´arquica. Mientras que en la multirresoluci´on uniforme los puntos de c´alculo de una rejilla son coincidentes (est´an autocontenidos) en los puntos de otra rejilla de mayor resoluci´on (figura 5.6), en la multirresoluci´on jer´arquica adaptativa los puntos de c´alculo son los puntos centrales de las celdillas, de manera que estos no ser´an coincidentes en los distintos niveles de resoluci´on, figura 5.13. El c´alculo del potencial y del error se realiza sobre las rejillas uniformes, como en la multirresoluci´on uniforme. Pero en un determinado nivel de resoluci´on, solo se calculan las celdillas que no est´an contenidas en una celda de la discretizaci´on jer´arquica correspondiente, celdas verdes en la figura 5.12. Las celdillas contenidas en otra de la discretizaci´on jer´arquica tendr´an su mismo valor de potencial, y no ser´a necesario calcular, celdillas rojas en la figura 5.12. El paso de una rejilla a la siguiente m´as fina Ωh ← Ω2h se realiza interpolando, ponderando con la distancia, en las celdillas que les corresponde subdivisi´on (figura 5.13) y manteniendo el mismo valor en las que deben permanecer sin refinamiento, celdillas rojas en la figura 5.12. En la figura 5.12 se ilustra los distintos niveles de discretizaci´on uniformes afectadas por los distintos niveles de discretizaci´on jer´arquica. I˜ niguez y Rosell en la publicaci´on [2] de 7.1, proponen un m´etodo de c´alculo en multirresoluci´on jer´arquica adaptativa.

5.3.1.

Interpolaci´ on

Al pasar de una discretizaci´on gruesa a la siguiente m´as fina, se obtienen los valores intermedios sobre la rejilla fina mediante una interpolaci´on (figura 5.13a). Esta interpolaci´on puede ser simplemente pasar el valor de la celda de menor resoluci´on a las de mayor resoluci´on, contenidas en ella (si no hay partici´on jer´arquica), o bien obtener los valores de las celdas realizando una interpolaci´on lineal ponderada, respecto de un par´ametro, de los valores de las celdas vecinas (si hay partici´on jer´arquica). Para obtener dicho par´ametro de ponderaci´on αi correspondiente a la celda ci , se considera la distancia eucl´ıdea d2 (c0 , ci ) = d21 + d22 + · · · + d2j + · · · + d2n (5.50) Donde dj = |qj (ci ) − qj (c0 )|,

(5.51)

As´ı, el par´ametro de ponderaci´on αi correspondiente a la celda ci , ser´a αi = D − d(ci , c0 ),

(5.52)

´ ´ CAP. 5. CALCULO NUMERICO

96

Figura 5.13: Ejemplos de interpolaci´on en una discretizaci´on jer´arquica adaptativa: (a) Interpolaci´on. (b) Mapeado. donde D vale, para n dimensiones:

D=

2n X

d(ci , c0 )

(5.53)

i=1

de manera que el peso de cada celda sea proporcional a su cercan´ıa al punto de c´alculo c0 . Con lo que se tiene: P2n α i ui , (5.54) u0 = Pi=1 2n i=1 αi

5.3.2.

Mapeado

Inversamente, al pasar de una discretizaci´on fina a la siguiente m´as gruesa, es necesario transferir o mapear los valores previamente calculados. Puesto que con la discretizaci´on en celdillas propuesta no coinciden los puntos, tambi´en ser´a necesario realizar una interpolaci´on lineal si existe subdivisi´on de celdas, en este caso la interpolaci´on ser´a no ponderada ya que los puntos involucrados ser´an equidistantes y por lo tanto no es necesario realizar c´alculo de distancias. Por ejemplo, el punto 4 de la figura 5.13b, como punto central de su celda, se computar´a como el valor medio de los cuatro puntos de sus subceldas, los cuales son equidistantes al mismo. En el caso de multirresoluci´on con rejillas uniformes los puntos de c´alculo se sit´ uan en los nodos de intersecci´on (figura 5.6) y son coincidentes en distintas resoluciones, por lo que el mapeado es directo sin necesidad de c´alculos. Sin embargo, el ahorro de puntos de c´alculo iterativo obtenido con la multirresoluci´on jer´arquica compensa el coste adicional del mapeado.

´ ´ CAP. 5. CALCULO NUMERICO

5.3.3.

97

Algoritmos

Los algoritmos necesarios para el c´alculo de la funci´on potencial, soluci´on de la ecuaci´on de Poison, se utilizar´an posteriormente en el cap´ıtulo 6 para la obtenci´on del planificador de trayectorias. 5.3.3.1.

Algoritmo de c´ alculo de la funci´ on potencial arm´ onica del C-espacio

Sobre una discretizaci´on jer´arquica adaptativa del C-espacio se realiza el c´alculo de la funci´on potencial utilizando el m´etodo de multirresoluci´on presentado en el subapartado (5.2.2). La funci´on base definida como - CalcularPotencialCespacioMV(nivelJM in, nivelJM ax) realizar´a un c´alculo de ida y vuelta en V entre un nivel de resoluci´on grueso definido por la variable nivelJM in y un nivel de resoluci´on fino definido por nivelJM ax (algoritmo 5.1). Estos dos niveles de resoluci´on estar´an comprendidos entre el m´ınimo y el m´aximo corrientes del proceso. Las funciones utilizadas por el algoritmo son: potencialRejilla(i, nivelJ): calcula en iter iteraciones el potencial de la rejilla de nivel jerarquico nivelJ, utilizando el m´etodo de relajaci´on S.O.R. mapearRejilla(nivelJ, nivelJ + 1): mapea el potencial de la rejilla de nivel (nivelJ) a la de nivel (nivelJ+1), siguiendo el procedimiento presentado en el subapartado (5.2.2.2). errorPotencialRejilla(i, nivelJM ax): calcular en i iteraciones el error sobre la rejilla (nivelJMax), utilizando las expresiones del residuo (5.42 y 5.43). interpolarRejilla(nivelJ, nivelJ−1): realiza una interpolaci´on para transferir los valores de la rejilla (nivelJ) a (nivelJ-1), tal como se muestra en el apartado (5.2.2.1). corregirRejilla(nivelJ − 1): corrige los valores interpolados sumando el error calculado (uk ← uk + ek ). De esta forma, la funci´on general - CalcularPotencialCespacioMJC (nivelJM in, nivelJM ax) utilizar´a la funci´on base, descrita en el algoritmo 5.1, para realizar el m´etodo recursivo denominado Multirresoluci´on Jer´arquica Completa (MJC) el cual efect´ ua un ciclo en V creciente (figura 5.8b), tal como se detalla en el apartado 5.2.

´ ´ CAP. 5. CALCULO NUMERICO

98

CalcularPotencialCespacioMV(nivelJMin,nivelJMax) begin for (nivelJ = nivelJM in to nivelJM ax − 1) do (P otencialRejilla) ← potencialRejilla(iter, nivelJ); (P otencialRejilla) ← mapearRejilla(nivelJ, nivelJ + 1); end; (errorP otencialRejilla) ← errorPotencialRejilla(iter, nivelJM ax); for (nivelJ = nivelJM ax to nivelJM in − 1) do (P otencialRejilla) ← interpolarRejilla(nivelJ, nivelJ − 1); (P otencialRejilla) ← corregirRejilla(nivelJ − 1); end; end; Algoritmo 5.1: C´alculo del potencial sobre una discretizaci´on jer´arquica adaptativa, utilizando el m´etodo de Multirresoluci´on en V.

CalcularPotencialCespacioMJC(nivelJMin,nivelJMax) begin (P otencialRejilla) ← inicializaP otencialAleatorioRejilla(nivelJM ax); for i = 1 to (nivelJM ax − nivelJM in) do (P otencialRejilla) ← CalcularP otencialCespacioM V (nivelJM ax − i, nivelJM ax); end; end; Algoritmo 5.2: C´alculo del potencial sobre una discretizaci´on jer´arquica adaptativa, utilizando el m´etodo de Full Multigrid (FMG). Donde: inicializaPotencialAleatorioRejilla(nivelJM ax) inicializa con valores aleatorios la variable vectorial PotencialRejilla en su m´axima resoluci´on. 5.3.3.2.

Algoritmo de obtenci´ on de la trayectoria

Una vez calculada la funci´on potencial, sobre la discretizaci´on jer´arquica, siguiendo el gradiente de la misma se obtiene el canal de trayectoria, el cual estar´a constituido por un conjunto de celdas concatenadas de diversos tama˜ nos que enlaza la celdilla de inicio (cini ) con la celdilla final (cf in ).

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.14: Funci´on potencial y canal de trayectoria.

Figura 5.15: Refinamiento del canal.

99

´ ´ CAP. 5. CALCULO NUMERICO

100

La funci´on TrayectoriaFina(T rayectoria), obtiene en refinamientos consecutivos la trayectoria en su m´axima resoluci´on (algoritmo 5.3), para ello se le pasa el canal a trav´es de la variable Trayectoria y devuelve en la misma varible la trayectoria deseada. Las funciones utilizadas por el algoritmo son: resolucion(T rayectoria): devuelve valor cierto cuando todas las celdillas que componen la variable Trayectoria son de la m´axima resoluci´on. particion(T rayectoria): devuelve las celdas de Trayectoria despu´es de realizar una partici´on 2n de aquellas que no est´en con la m´axima resoluci´on. calculoPotencial(Canal): realiza y devuelve el c´alculo del potencial en el dominio definido por las celdas de la variable vectorial Canal, considerando contorno la parte exterior del mismo. canalTrayectoria(Canal, bini , bf in ): devuelve el nuevo canal de trayectoria obtenido en el dominio formado por Canal, con la celda inicial y final indicadas. Trayectoria(Canal, bini , bf in ) begin while ¬resolucion(T rayectoria) do begin (Canal) ← particion(T rayectoria); (Canal) ← calculoPotencial(Canal); (T rayectoria) ← canalTrayectoria(Canal, bini , bf in ); end; return T rayectoria end; Algoritmo 5.3: Algoritmo para la obtenci´on de la trayectoria a partir del canal.

5.4.

Obtenci´ on de resultados

La inexistencia de aplicaciones que permitan implementar los algoritmos presentados en este cap´ıtulo, ha hecho necesario desarrollar una mediante la cual se puede obtener resultados num´ericos y gr´aficos en distintos entornos. La herramienta de programaci´on utilizada ha sido Visual C++ de Microsoft. • En la figura 5.14 se muestra, en 2D, un ejemplo de entorno discretizado jer´arquicamente, con la funci´on potencial calculada modelando los obst´aculos como contorno, para unas condiciones de contorno de Neuman y suponiendo continuidad c´ıclica de la variable q1 del espacio de configuraciones.

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.16: Aplicaciones distintos escenarios.

101

102

´ ´ CAP. 5. CALCULO NUMERICO

• En la figura 5.15 se muestra un ejemplo en 2D, de los pasos seguidos (de izquierda a derecha y de arriba hacia abajo) para la obtenci´on de una trayectoria en m´axima resoluci´on, partiendo del canal de trayectoria. • En la figura 5.16 se muestra, tambi´en en 2D, la obtenci´on de trayectorias en distintos entornos con grados de dificultad extremos (excepto el caso a), con lo que se puede corroborar el alto grado de eficiencia del m´etodo en lo que se refiere a la obtenci´on de trayectorias. En los cuatro casos planteados (a, b, c y d) se modelan los obst´aculos como contorno, el punto inicial como valor elevado y el punto final como valor bajo y utilizan condiciones de contorno de Neuman. Se muestra primero la figura del entorno (completamente discretizado) con los dos puntos extremos de la trayectoria requerida, la siguiente figura incorpora el canal de trayectoria encontrado y la tercera figura muestra la trayectoria obtenida por refinamiento del canal. En este ejemplo se constata la conveniencia de utilizar condiciones de contorno de Neuman en los pasillos largos y estrechos, tal como se analiza en el apartado (2.3.2). • En la figura 5.17 se muestra, gr´aficamente, los resultados computacionales (n´ umero de iteraciones) obtenidos para tres tipos de entornos planteados: (a) entorno muy fragmentado, (b) medianamente fragmentado y (c) poco fragmentado con pasillos extrechos. El metodo iterativo utilizado es el S.O.R. estudiado en el apartado 5.1.4.3, donde con un entorno regular simple se obtuvo una (ω) ´optima. Sin embargo en este ejemplo se estudia su influencia para distintos entornos complicados, por lo que se hace variar entre cero y uno a intervalos de una cent´esima. Las condiciones de contorno planteadas ser´an mixtas ∂u u(r) = αuD + (1 − α) (r), ∀r ∈ ∂Ω (5.55) ∂n donde el par´ametro α se hace variar tambi´en entre cero y uno, a intervalos de una cent´esima. As´ı, un punto de la superficie indica el n´ umero de iteraciones realizadas para encontrar una trayectoria, con un determinado valor de α y de ω. El n´ umero de puntos de la superficie (104 puntos) representa el n´ umero de realizaciones o trayectorias encontradas. Entonces, se observa que existe un valor ´optimo del par´ametro ω (situado en la parte m´as baja del valle de las figuras 3D), el cual var´ıa ligeramente con respecto a α y tambi´en con respecto al tipo de entorno. Tambi´en se observa, con respecto a las condiciones de contorno, una situaci´on ´optima para (α = 1) correspondiente a las condiciones de contorno de Neuman, conforme α disminuye aumentan el n´ umero de iteraciones hasta alcanzar un m´aximo, disminuyendo despu´es con una pendiente que depende del tipo de entorno, pero que no llega hasta el m´ınimo de Neuman. • En la figura 5.18 se muestra la discretizaci´on jer´arquica de un entorno que modeliza los obst´aculos como funciones real positivas con forma de coronas y en torno al punto destino, fijado a un valor bajo, se sit´ ua una corona atractiva como una funci´on real negativa, adem´as se consideran condiciones de contorno en el infinito. En el apartado 3.7.3 se propone este enfoque

´ ´ CAP. 5. CALCULO NUMERICO

Figura 5.17: Iteraciones distintos escenarios.

103

´ ´ CAP. 5. CALCULO NUMERICO

104

Figura 5.18: Planificador mediante anillos. con el que se obtiene una distribuci´on del gradiente que no presenta zonas excesivamente planas.

5.5.

Aportaci´ on

En este cap´ıtulo se han estudiado los m´etodos num´ericos para el c´alculo de las funciones arm´onicas, subarm´onicas y superarm´onicas: - Se han estudiado los distintos m´etodos de relajaci´on, as´ı como la relaci´on existente entre la velocidad de convergencia de los mismo y el grado de resoluci´on. - En este cap´ıtulo tambi´en se ha propuesto un m´etodo num´erico novedoso para el c´alculo de las funciones arm´onicas, subarm´onicas y superarm´onicas, basado en m´etodos de multirresoluci´on aplicados a una discretizaci´on jer´arquica adaptativa. Con este m´etodo num´erico se obtiene una reduci´on de ciclos de iteraci´on al mismo tiempo que una reduci´on de puntos de c´alculo, todo ello conservando la m´axima resoluci´on en los contornos; condici´on necesaria para obtener funciones potenciales que sirvan para generar trayectorias pr´oximas a los obst´aculos, tal como ser´ıa el caso en pasillos estrechos.

Cap´ıtulo 6 Uso de las funciones arm´ onicas en espacios de configuraciones explorados mediante t´ ecnicas de muestreo Con la discretizaci´on jer´arquica adaptativa, introducida en el apartado 5.2, se consigue una reducci´on del n´ umero de celdillas de c´alculo respecto de la discretizaci´on uniforme. Esta reducci´on es tanto menor cuanto mayor es el contorno de obst´aculos (figura 6.1a, b y c) y, el n´ umero de celdillas totales crecer´a exponencialmente con el n´ umero de grados de libertad del C-espacio. Por este motivo, es costoso computacionalmente el calcular de forma exacta o aproximada (con cierta resoluci´on) todo el espacio de configuraciones; mientras que los m´etodos de muestreo son eficientes ya que no modelan el C-espacio, solo captan con curvas unidimensionales la conectividad del espacio libre del C-espacio. En esta tesis se presenta un nuevo planificador basado en la combinaci´on de m´etodos de muestreo con los m´etodos basados en funciones arm´onicas. Esta combinaci´on se basa en un entrelazado entre el c´alculo de funciones arm´onicas sobre un espacio discretizado jer´arquicamente y la exploraci´on de

Figura 6.1: Discretizaci´on jer´arquica completa del C-espacio en 2D: Efecto del crecimiento del contorno sobre la discretizaci´on. 105

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

106

dicho espacio mediante muestreo. Para ello se utiliza un tipo de muestreo que, en lugar de muestrear si una configuraci´on es libre o de colisi´on, se muestrea si una celdilla de discretizaci´on es libre, de obst´aculo o de contorno. As´ı, en el apartado 6.1 se presenta un planificador que utiliza una funci´on detectora de colisi´on que tambi´en devuelve la distancia de colisi´on, mientras que en apartado (6.4) se presenta una alternativa de planificador, que utiliza una funci´on detectora de colisi´on que no devuelve la distancia de colisi´on.

6.1.

Planificador PHM b´ asico

El m´etodo de planificaci´on propuesto, denominado “Probabilistic Harmonic Funtion Method”(PHM), resuelve la trayectoria utilizando una funci´on potencial arm´onica calculada sobre una discretizaci´on jer´arquica adaptativa del C-espacio; el cual es explorado de manera aleatoria con el nuevo tipo de muestreo de celdas, propuesto en este apartado, con el que se obtiene un conocimiento parcial del entorno en zonas y en resoluci´on. Los obst´aculos forman parte del contorno, sobre el que se aplica condiciones de contorno Dirichlet, y en el punto destino se sit´ ua la funci´on delta. Inicialmente el espacio de configuraciones es desconocido; u ´nicamente un punto de situaci´on y otro de destino (figura 6.2a), ambos de m´axima resoluci´on, se establecen determinando una discretizaci´on inicial que reporta una serie de celdas de distinto tama˜ no y contenido tambi´en desconocido. El muestreo aleatorio de estas celdas definir´a la situaci´on y forma de los obst´aculos progresivamente. Suponiendo que se muestrean las celdas c1 , c2 y c3 (figura 6.2b), estas se identifican como libre, obst´aculo y de contorno respectivamente; c3 se subdivide generando 2n celdas que se a˜ naden al conjunto de celdas de contenido desconocido (figura 6.2c). De esta forma, en la discretizaci´on del C-espacio, se establecen tres tipos o conjuntos de celdas: Celdas de contenido desconocido (grises) G = {cG1 , cG2 , . . . , cGNG },

(6.1)

Celdas pertenecientes a los obst´aculos (negras) B = {cB1 , cB2 , . . . , cBNB }

(6.2)

Celdas pertenecientes al espacio libre (blancas) W = {cW1 , cW2 , . . . , cWNW }

(6.3)

Los conjuntos de celdas G, B y W se definir´an como variables vectorizadas, que estar´an formadas por los c´odigos o etiquetas de las celdas correspondientes, pudiendo ser estas de m´axima resoluci´on o no, dependiendo del nivel de discretizaci´on.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

107

Figura 6.2: Muestreo y clasificaci´on de celdas del C-espacio en 2D. El muestreo aleatorio se realiza sobre el conjunto de celdas grises (G) cuyo n´ umero, partiendo de un n´ umero (NG ) reducido (figura 6.2a), crece inicialmente y posteriormente se va reduciendo, mientras los C-obst´aculos se van definiendo progresivamente. Sobre este escenario discretizado jer´arquicamente, se aplican las condiciones Sorteo(G) begin F ← FuncionIntegral (G); r ← random (0, 1); {cGk ← G} tal que {F (cGk ) ≤ r}; return cGk end; Algoritmo 6.1: Algoritmo del sorteo para el muestreo aleatorio de una celda perteneciante a G.

de Dirichlet y se calcula una funci´on potencial V (apartado 5.2.4), que se utiliza para buscar un canal de trayectoria; lo cual puede ocurrir con un conocimiento parcial del entorno y una reducci´on del n´ umero de puntos de c´alculo. Con este enfoque se puede obtener un planificador completo, probabil´ısticamente y en resoluci´on. Dicho planificador se desarrolla en los subapartados siguientes, y es presentado por I˜ niguez y Rosell en la publicaci´on [3] del apartado 7.1.

6.1.1.

Selecci´ on aleatoria de una celda ponderada por su tama˜ no: sorteo para la exploraci´ on de celdas

Sea cGk una celda perteneciente al conjunto G, la probabilidad de selecci´on de cGk vendr´a determinada por el peso ω(cGk ), que tendr´a en cuenta el tama˜ no de dicha celda. Suponiendo que cGk es una celda de nivel de partici´on mk ,

108

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

Figura 6.3: Muestreo aleatorio utilizando la funci´on integral F (ck ). entonces el peso ω(cGk ) se expresa como: ω(cGk ) =

VcGk VT

=

1 2mk n

,

(6.4)

n y VT es el voludonde VcGk es el volumen de la celdilla cGk en la rejilla Gm n men total de todas las cedillas perteneciente a Gm . Es decir, la probabilidad de selecci´on de una celda se incrementar´a con el tama˜ no de la misma; de esta forma se producir´a una r´apida caracterizaci´on de los contornos de los obst´aculos, ya que la incertidumbre de las celdas grandes parcialmente libres ser´a despejada antes. El sorteo efectuado en la operaci´on de muestreo se realiza utilizando una funci´on integral F (cGk ) definida sobre el conjunto de celdas G

F (cGk ) =

k X i=1

ω(cGi ) ∀cGi ∈ G y

1 ≤ k ≤ NG

(6.5)

As´ı, el valor m´aximo de esta funci´on integral se obtendr´a para la u ´ltima celda de G (cGN G ), y valdr´a: Fmax = F (cGN G ) = 1

(6.6)

El procedimiento de muestreo aleatorio o sorteo se implementa de la manera siguiente: Se obtiene un n´ umero aleatorio (r), comprendido entre cero y uno, utilizando la funci´on integral F (cGk ) se extrae la celdilla correspondiente (cGk en el ejemplo de la figura 6.3 y algoritmo 6.1).

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

109

Figura 6.4: Clasificaci´on de las celdas en el espacio de configuraciones.

6.1.2.

Verificaci´ on y clasificaci´ on de celdas

La celda cGk , extraida por sorteo (con ponderaci´on del tama˜ no) de G, se verifica y se identifica como celda libre, obst´aculo o contorno. La configuraci´on del robot correspondiente a su centro se pasa a la funci´on detectora de colisi´on, la cual reportar´a la distancia de colisi´on (dk ) sobre el espacio de configuraciones (esta distancia ser´a positiva si la configuraci´on corresponde al espacio libre y negativa si corresponde al espacio obst´aculo). Entonces, de acuerdo con la figura 6.4, la celda se clasifica como: Celda libre: Celda inscrita en la circunferencia de radio dk y dk > 0 (celda c1 ) Celda obst´aculo: Celda inscrita en la circunferencia de radio dk y dk < 0 (celda c2 ) Celda contorno: Celda no inscrita en la circunferencia de radio dk (celda c3 ) Las celdas clasificadas como libre y obst´aculo se clasifican como celdas W y B respectivamente, mientras que la celda contorno se subdivide y las 2n celdas resultantes se clasifican como celdas de contenido desconocido G. El apartado 6.2 describe este procedimiento, usando las siguientes funciones: DistanciaColision(cG ): devuelve la distancia m´ınima (en el C-espacio) del centro de la celda cG al contorno del obst´aculo m´as proximo. Particion(cG ): devuelve la partici´on de la celda cG .

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

110

VerificaClasifica (cG ) begin d ← DistanciaColision(cG ); if (|d| ≥ D) then if (d > 0) then (W ← W ∪ cG ); else (B ← B ∪ cG ); else { (cG1 , cG1 , . . . , cGH ) ← Particion (cG ); G ← G ∪ (cG1 , cG1 , . . . , cGH ); } end;

Algoritmo 6.2: Algoritmo para la verificaci´on y clasificaci´on de celdas. D es la distancia diagonal de la celda a verificar.

6.1.3.

Algoritmo del planificador de trayectorias. Exploraci´ on pasiva de trayectorias

Combinando el m´etodo anterior de exploraci´on de celdas, el c´alculo de funciones potenciales arm´onicas sobre una discretizaci´on jer´arquica y la exploraci´on pasiva de trayectorias (consistente en calcular trayectorias sobre zonas parcialmente desconocidas del C-espacio y evaluaci´on posterior) se dise˜ na el generador de trayectorias (algoritmo 6.3): • Inicialmente, se establece la posici´on de los puntos inicial y final definidos por sendas celdas de m´axima partici´on, cini y cf in . Con el resto del C-espacio se realiza una discretizaqci´on jer´arquica (figura 6.5a), clasificando todas las celdas como de contenido desconocido (grises), excepto cini y cf in . • Despu´es de la exploraci´on de NE celdas (figura 6.5b), se aplican las condiciones de contorno de Dirichlet en los obst´aculos (celdas negras) y se calcula la funci´on potencial (VW G ) sobre la discretizaci´on actual, considerando las celdas grises como libres, y se realiza la b´ usqueda de un canal trayectoria (figura 6.5c). • Las celdas grises pertenecientes a este canal de trayectoria, ser´an

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

111

PlanificadorTrayectoria(cini ,cf in ) begin (cini , cf in ) ← MarcaInicioFin; G ← DiscreInicial (cini , cf in ); B ← ∅; W ← ∅; T rayectoria ← ∅; i ← 0; do { for j = 1 to NE (cG ) ← Sorteo(G); VerificaClasifica(cG ); end for VW G ← CalculaPotencial(G, W, B, cf in ); Canal ← CanalTrayectoria(VW G ); Libre ← CompruebaClasifica(Canal); if (Libre == T RU E) T rayectoria ← Trayectoria(Canal); return T rayectoria; end if ++i; } while ((Libre == F ALSE) AN D (i < NI )); return T rayectoria; end; Algoritmo 6.3: Algoritmo planificador de trayectorias, PHM b´asico.

comprobadas y clasificadas. Si el canal de trayectoria queda interrumpido, por celdas subdivididas y clasificadas como grises, se reinicia el procedimiento. Repitiendo el proceso hasta que se encuentre un canal de trayectoria libre o se hayan realizado un n´ umero de iteraciones dado (NI ). • Por u ´ltimo, el canal de trayectoria encontrado se refina hasta obtener la trayectoria con la m´axima resoluci´on (figura 6.5d). Para ello se sigue el m´etodo presentado en el subapartado 5.2.3.2. Para realizar estas tareas el algoritmo cuenta con las siguientes funciones: CalculaPotencial(G, W, B, cf in ): calcula la funci´on potencial considerando las celdas grises junto con las celdas libres. CanalTrayectoria(VW G ): utilizando una funci´on de navegaci´on (por seguimiento del gradiente del potencial) intenta encontrar un canal de trayectoria: si lo encuentra devuelve los c´odigos de las celdas que lo componen en la variable vector Canal y si no devuelve (Canal = ∅).

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

112

Figura 6.5: Secuencia ejemplo del m´etodo PHM b´asico. CompruebaClasifica(Canal): cada una de las celdas grises que forman parte del canal son verificadas: si son todas libres devuelve (Libre = T RU E) y si no devuelve (Libre = F ALSE). DiscreInicial (cini , cf in ): partiendo del conocimiento de la posici´on de las celdas inicial y final, con la m´axima resoluci´on, obtiene una discretizaci´on jer´arquica inicial del todo el espacio, en celdas grises de contenido desconocido (figura 6.5a) Trayectoria(Canal): partiendo del Canal de trayectoria obtiene la trayectoria compuesta por celdillas de la m´axima resoluci´on (figura 6.5d)

6.1.4.

Estudio de la completitud del m´ etodo propuesto

Desde un punto de vista din´amico pueden verse las celdas contorno de un determinado nivel de partici´on como poblaciones de generaciones sucesivas (figura 6.6), de forma que la poblaci´on de la generaci´on correspondiente al nivel m dar´a lugar a la poblaci´on de la generaci´on siguiente (m + 1). El simil ser´ıa: una celda-individuo seleccionado al azar, el cual pertenece a una generaci´on determinada m, genera 2n hijos pertenecientes a la generaci´on (m + 1) y muere. En un instante determinado coexisten distintas generaciones pr´oximas, creciendo en poblaci´on la generaci´on m´as joven mientras se van extinguiendo las m´as viejas. Al mismo tiempo, se observa un crecimiento demogr´afico global hasta alcanzar un m´aximo, descendiendo despu´es hasta la extinci´on total de la especie ya que la generaci´on de celdas de m´axima resoluci´on mueren sin descendencia.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

113

Figura 6.6: Ejemplo de evoluci´on de la poblaci´on de celdas en G. Este estudio se desarrolla en los subapartados siguientes, y es presentado por I˜ niguez y Rosell en la publicaci´on [3] del apartado 7.1. 6.1.4.1.

Modelo de poblaci´ on de celdas en G

La probabilidad de que una determinada m-celda sea explorada es obtenida estudiando la evoluci´on de las poblaciones de las celdas parcialmente libres en sus distintos niveles jer´arquicos. En primer lugar se obtiene el modelo de crecimiento de poblaciones para una determinada tasa de nacimientos y muertes, para ello se define la nomenclatura siguiente: • xm,k : es la poblaci´on o n´ umero de celdas de nivel m en G despu´es de k celdas exploradas. • Pm : es la probabilidad de que una m-celda sea parcialmente libre. • Pxm,k : es la probabilidad de que en el paso k el algoritmo seleccione la poblaci´on xm,k . • Pm,k : es la probabilidad de que en el paso k el algoritmo seleccione una determinada celda perteneciente a la poblaci´on xm,k . • Pm,1..k : es la probabilidad de que en k pasos el algoritmo seleccione una determinada celda perteneciente a la poblaci´on xm,k . Cuando una m-celda de G es seleccionada y evaluada como parcialmente libre, entonces se subdivide en 2n (m + 1)-celdas hijas y se retornan a G, elimin´andose la m-celda progenitora. Por lo tanto se cumplen las siguientes consideraciones: • xm,k : decrementa en uno si en el paso (k − 1) una m-celda ha sido seleccionada.

114

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

Figura 6.7: Modelo de poblaci´on de celdas en G.

• xm,k : incrementa en 2n si en paso (k − 1) una (m-1)-celda ha sido seleccionada y clasificada como parcialmente libre. • xm,k : permanece constante si en el paso (k − 1) no han sido seleccionadas una m-celda ni una (m-1)-celda, o s´ı ha sido seleccionada una (m-1)-celda y no era parcialmente libre . Por lo tanto, la poblaci´on o n´ umero de celdas parcialmente libres de nivel jer´arquico (generaci´on) m vendr´a dada por el siguiente modelo probabil´ıstico: xm,k = (xm,k−1 − 1) · Pxm,k +

(xm,k−1 + 2n ) · Pxm−1,k · Pm−1 +

xm,k−1 · (1 − Pxm,k − Pxm−1,k )+

xm,k−1 · Pxm−1,k · (1 − Pm−1 ) (6.7)

En la figura 6.6 se muestra un ejemplo de evoluci´on de poblaciones y en la figura 6.7 se muestra la estimaci´on de poblaciones mediante el modelo probabil´ıstico propuesto. Considerando un espacio de configuraciones de n dimensiones con una discretizaci´on jer´arquica de M niveles, el n´ umero de celdas existentes en un nivel m tal que m ≤ M es de 2mn m-celdas. Si Nm es el n´ umero de m-celdas que contienen parte del contorno, entonces la probabilidad de que una m-celda sea parcialmente libre es: Nm (6.8) Pm = nm 2 Este valor ser´a estimado muestreando un n´ umero sm de m-celdas y verificando el n´ umero de celdas gm parcialmente libres encontradas, entonces Pm se estima como: gm (6.9) Pm ≃ sm Las M-celdas que son parcialmente libres son consideradas como celdas obst´aculo, puesto que son de m´axima resoluci´on, por lo tanto se cumple: PM = 0

(6.10)

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

115

Figura 6.8: Probabilidad Pxm,k de seleccionar la poblaci´on xm,k para m = 3, 4, 5 y 6. 6.1.4.2.

Probabilidad de selecci´ on condicionada por el tama˜ no de las celdas

En un espacio discretizado jer´arquicamente el tama˜ no de las celdas es heterog´eneo y depende del nivel jer´arquico de partici´on. Para una C-espacio de n dimensiones, si se considera unitario el tama˜ no del hipervolumen de una M -celda, el tama˜ no de una celda de nivel de partici´on m viene dado por el hipervolumen Vm = 2n(M −m) (6.11) Por lo tanto el muestreo aleatorio debe ser hecho considerando el tama˜ no de las celdas. Primero se selecciona la poblaci´on de un determinado nivel jer´arquico y despu´es se escoge una celda de esta generaci´on para ser explorada: a) Sea Vm,k el volumen de todas las m-celdas de G en el paso k. La probabilidad Pxm,k depende del volumen Vm,k con respecto al total del volumen de celdas en G: Vm,k 2n(M −m) · xm,k 2−nm · xm,k Pxm,k = PM = PM = PM n(M −j) · x −nj · x j,k j,k j=1 Vj,k j=1 2 j=1 2

(6.12)

En la figura 6.8 se muestra un ejemplo de la evoluci´on real de Pxm,k y en la figura 6.9 su estimaci´on, utilizando el modelo de poblaci´on (ecuaci´on 6.7) b) La probabilidad de seleccionar una m-celda de xm,k es: Pm,k =

1 xm,k

(6.13)

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

116

Figura 6.9: Estimaci´on de la probabilidad Pxm,k para m = 4, 5 y 6. El producto de las dos anteriores da la probabilidad de que en el paso k el algoritmo seleccione una m-celda de G. De aqu´ı, la probabilidad de explorar una determinada m-celda de G en k tiradas es (figura 6.10): Pm,1..k = 1 −

6.1.4.3.

k Y j=1

1 − Pxm,j Pm,j




(6.14)

Probabilidad de fallo

Asumiendo que existe una ruta entre dos configuraciones ci y cf formada por una cadena de celdas cuyos tama˜ nos ser´an iguales o superiores al de una M-celda, el objetivo es determinar el n´ umero m´ınimo k de iteraciones del algoritmo tal que la probabilidad de fallo del planificador permanezca dentro de un l´ımite dado. P umero total de Sea Dm el n´ umero de m-celdas del canal, y D = Lm=1 Dm el n´ celdas del canal. Entonces la probabilidad de fallo, despu´es de k iteracciones, ser´a: M Y (6.15) (Pm,1..k )Dm P (F allo) = 1 − m=1

Puesto que la probabilidad de seleccionar una m-celda del canal disminuye con el tama˜ no de la misma, se considera el peor caso en el que el canal estar´a formado por DM M-celdas. Entonces, la probabilidad de fallo ser´a: P (F allo) ≤ 1 − (PM,1..k )DM

(6.16)

Para una determinada cota de fallo (P (F allo) ≤ δ), se deber´a cumplir PM,1..k ≥ (1 − δ)(1/DM )

(6.17)

El n´ umero de iteraciones k del algoritmo se obtiene aplicando iterativamente la ecuaci´on (6.14) para (m = M ) hasta que se cumpla la desigualdad anterior.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

117

Figura 6.10: Probabilidad de seleccionar una m-celda en k intentos (ecuaci´on (6.14)). La curva correspondiente al nivel M es utilizada para obtener el n´ umero de iteraciones necesarias.

Figura 6.11: Ejemplo de C-espacio para dos n´ umeros de muestras: a) k = 825 (exploraci´on exhaustiva) y b) k=80. 6.1.4.4.

Validaci´ on

Los resultados experimentales se han obtenido sobre un espacio de configuraciones de dos dimensiones, con dos obst´aculos y una discretizaci´on no regular, jer´arquica adaptativa con (M = 6) niveles (figura 6.11). En la figura 6.11a se han realizado el n´ umero de muestras necesarias para efectuar una exploraci´on exhaustiva, mientras que en 6.11b se ha realizado una exploraci´on parcial obteni´endose celdas libres, obst´aculo y parcialmente libres. Suponiendo una cota de fallo (δ = 0,1) y dos configuraciones ci y cf para ser conectadas a trav´es de un canal compuesto por (D = 20) celdas, la probabilidad de que una m-celda sea parcialmente libre es estimada como: m Pm =

gm sm

1 1

2 0.66

3 0.69

4 0.59

5 0.45

6 0

Usando estos valores, la figura 6.7 muestra la poblaci´on estimada de m-celdas parcialmente libres computadas utilizando la ecuaci´on (6.7), la figura 6.9

118

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

Figura 6.12: Canal de trayectoria con exclusi´on de celdas. muestra la evoluci´on de las probabilidades estimadas utilizando la ecuaci´on (6.12). El n´ umero de iteraciones debe ser tal que se cumpla: Pm,1..k ≥ (1 − δ)(1/D) = 0,91/20 = 0,995

(6.18)

Aplicando iterativamente la ecuaci´on (6.14), esta condici´on se satisface para k ≥ 655 (figura 6.10). Utilizando este valor de k, el algoritmo 6.3 ( PlanificadorTrayectoria) se ha utilizado para explorar el C-espacio. El resultado se muestra en la figura 6.11b donde la funci´on arm´onica se ha calculado sobre libres y se ha utilizado para encontrar la trayectoria entre ci y cf .

6.2.

Alternativas de mejora del m´ etodo b´ asico

Sobre el esquema b´asico del planificador propuesto en el apartado (6.1) caben m´etodos de mejora basados en reglas heur´ısticas. La aplicaci´on de estos m´etodos permite la realizaci´on de un planificador m´as eficiente en tiempo de ejecuci´on. En los subapartados siguientes se proponen varias posibilidas de mejora (algunas ya estudiadas) y se aplican en el apartado (6.3), en el dise˜ no del algoritmo de un planificador avanzado. En las publicaciones [4 y 5] del apartado 7.1, Rosell e I˜ niguez proponen alternativas de mejora del planificador PHM b´asico.

6.2.1.

Exclusi´ on de celdas del sorteo de exploraci´ on

Considerando los obst´aculos como contorno y utilizando condiciones de contorno de Dirichlet, el valor de VW G se incrementa monot´onicamente desde el punto final hasta los C-obst´aculos. Solamente aquellas celdas con un valor de VW G m´as bajo que el de la celda inicial deben ser consideradas para el muestreo, ya que el canal de trayectoria es encontrado siguiendo el gradiente de la funci´on VW G . De esta forma, aquellas celdas de G con valores de potencial superiores al de la celda inicial pueden excluirse del sorteo, reduci´endose

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

119

el tiempo de exploraci´on. Es decir, aquellas celdas que cumplen: cGi ∈ G | VW G (cGi ) > VW G (cini ),

(6.19)

son transferidas al conjunto A = {cA1 , cA2 , . . . , cAN A },

(6.20)

quedando excluidas del sorteo de muestreo, consider´andose como celdas obst´aculo. De esta manera se acelera el proceso; a no ser que estas nuevas celdas bloqueen un posible camino hacia el objetivo y por tanto el algoritmo no encuentre un camino libre, en cuyo caso todas las celdas del conjunto A se devuelven al conjunto G, continuando la b´ usqueda. En la figura 6.12 se muestra un ejemplo (para 2D) donde se ha encontrado una trayectoria despu´es de haber excluido un n´ umero de celdas (azules) que cumpl´ıan la condici´on (6.19).

6.2.2.

Exploraci´ on pasiva de trayectorias: Puntos de Apoyo y Funci´ on Potencial Auxiliar

En el apartado (6.1.3) se present´o el m´etodo para la exploraci´on de trayectorias: considerando las celdas grises como libres se obtiene una trayectoria, y aquellas celdas grises pertenecientes a dicha trayectoria ser´an comprobadas a posteriori. Si la trayectoria queda interrumpida por celdas subdivididas y clasificadas como grises, las celdas extremos de los trozos de trayectoria resultantes ser´an consideradas celdas de apoyo para la exploraci´on subsiguiente, clasific´andose dentro del conjunto Y Y = {cY1 , cY2 , . . . , cYN Y },

(6.21)

que incluir´a tambi´en las celdas inicial y final. Estas celdas de apoyo ser´an prefijadas a valores bajos de una nueva funci´on potencial VW GY = VW GY (ck ) | ck ∈ (W ∪ G), (6.22) la cual se utilizar´a para modular las probabilidades de selecci´on en el sorteo de exploraci´on de celdas y que ser´a soluci´on de la ecuaci´on de Laplace; considerando los obst´aculos como contorno y aplicando condiciones de contorno de Dirichlet, con valor alto para el contorno (Vcontorno ) y bajo para los puntos de apoyo (Vapoyo ). Si la probabilidad de selecci´on de una celda es inversamente proporcional al valor de la funci´on potencial VW GY , entonces los entornos de los puntos de apoyo ser´an explorados antes, conduci´endose la exploraci´on hacia aquellas zonas m´as prometedoras de contener caminos libres. Con este planteamiento, la funci´on integral para el sorteo de celdas depender´a del peso ω(ck ) = ω1 (ck ) · ω2 (ck ), (6.23)

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

120

donde ω1 (ck ) es el peso correspondiente al tama˜ no (ecuaci´on (6.4)) y ω2 (ck ) =

Vcontorno − VW GY (ck ) , Vcontorno − Vapoyo

(6.24)

que tendr´a un valor comprendido entre cero (para las celdas de contorno) y uno (para las celdas de apoyo). La funci´on integral para el sorteo de celdas, se obtiene como en la ecuaci´on (6.5):

F (ck ) =

k X i=1

ω(ci ) ∀ci ∈ G y

1 ≤ k ≤ NG

(6.25)

Por u ´ltimo, con la finalidad de aprovechar la recursividad en el c´alculo de las funciones arm´onicas, se memorizar´an en variables distintas las funciones potenciales calculadas para la obtenci´on de trayectorias (VW G ) y las calculadas para obtener probabilidades (VW GY ).

6.2.3.

Resoluci´ on m´ axima ajustable

Teniendo en cuenta que, en la discretizaci´on jer´arquica adaptativa del C-espacio, es posible prefijar la m´axima resoluci´on alcanzable (figura 5.12) y teniendo en cuenta que las posibles trayectorias no siempre necesitar´an de una resoluci´on m´axima lo m´as fina posible; entonces resultar´a eficiente, respecto a la reducci´on de puntos de c´alculo, el ajuste pasivo del par´ametro identificativo de la m´axima resoluci´on (M ). Con la finalidad de reducir al m´aximo el n´ umero de puntos de c´alculo, la m´axima resoluci´on alcanzable en el proceso de discretizaci´on jer´arquica se aplica con un criterio basado en reglas: La m´axima resoluci´on posible se aplicar´a siempre a los puntos inicial y destino. La m´axima resoluci´on aplicable a los contornos ser´a variable por fases de b´ usqueda, comenzando por una resoluci´on m´as gruesa y afinando seg´ un va abanzando el proceso sin obtenci´on de resultado.

6.2.4.

Aplicaci´ on de condiciones de contorno de Neuman en los obst´ aculos

La aplicaci´on de condiciones de contorno de Dirichlet en los l´ımites del C-espacio y de Neuman en los obst´aculos (figuras 5.17 y 5.18) consigue que los pasillos formados entre ellos no se cieguen, conservando un gradiente longitudinal m´ınimo que permita el paso de trayectorias.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

6.3.

121

Algoritmo del planificador PHM avanzado

Aplicando las mejoras del apartado anterior al planificador PHM b´asico se obtiene el planificador PHM avanzado (algoritmo 6.5), donde la b´ usqueda de trayectoria es conducida por zonas m´as prometedoras, y donde se produce una reduci´on de puntos de c´alculo. Mientras que los puntos b´asicos del algorito son conservados, se destacan en negrita los puntos nuevos de mejora: • Inicialmente, se establece la posici´on de los puntos inicial y final definidos por sendas celdas de m´axima partici´on, cini y cf in . Con el resto del Cespacio se realiza una discretizaci´on jer´arquica, clasificando las celdas como de contenido desconocido. Adem´ as se establece un ´ındice de m´ axima resoluci´ on (M ) que, en principio, no ser´ a el m´ aximo posible. • El conjunto de celdas de apoyo estar´ a formado inicialmente por las celdas cini y cf in . • Las celdas de apoyo ser´ an prefijadas a valores bajos de potencial y se calcular´ a la funci´ on potencial auxiliar (VW GY ), la cual se utilizar´ a para modular las probabilidades de selecci´ on en la pr´ oxima exploraci´ on. • Despu´es de la exploraci´on de NE celdas, se aplican las condiciones de contorno de Dirichlet en la frontera del C-espacio, Neuman en los obst´aculos y se calcula la funci´on potencial (VW G ) sobre la discretizaci´on actual, consider´andose las celdas grises como libres, y realiz´andose la b´ usqueda de un canal trayectoria. • Las celdas grises pertenecientes a este canal de trayectoria, ser´an comprobadas y clasificadas. Si el canal de trayectoria queda interrumpido por celdas subdivididas y clasificadas como grises, las celdas extremo de los trozos resultantes ser´ an consideradas celdas de apoyo. • Utilizando la funci´ on potencial VW G y teniendo en cuenta la ecuaci´ on (6.19), se obtiene el conjunto (A) de celdas que se excluir´ an del pr´ oximo sorteo. • Si no se obtiene el canal de trayetoria y el conjunto A es no vac´ıo; entonces se restituyen los elementos de dicho conjunto al de celdas grises. • Despu´ es se incrementa el ´ındice de m´ axima partici´ on (M ), si se ha llegado al m´ aximo (MM ) se aplica un anillo de atracci´ on en torno del punto destino y se vuelve al principio, repiti´ endose el proceso hasta que se encuentre un canal de trayectoria libre o se

122

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

hayan realizado un n´ umero de iteraciones dado (NI ). El procedimiento detallado se muestra en el algoritmo (6.5), donde las nuevas funciones incluidas en el algoritmo del planificador PHM avanzado son: CalculaPotencial(G, W, B, Y ): calcula la funci´on potencial auxiliar (VGW Y ) considerando las celdas grises junto con las celdas libres y las celdas de apoyo incluidas en el conjunto Y . ObtencionCeldasA(VW G ): obtiene las celdas que ser´an incluidas en el conjunto A y que ser´an excluidas del pr´oximo sorteo. PlanificadorTrayectoria(cini ,cf in ) begin (cini , cf in ) ← MarcaInicioFin; G ← DiscreInicial (cini , cf in ); B ← ∅; W ← ∅; A ← ∅; P ← ∅; Y ← {cini , cf in }; T rayectoria ← ∅; i ← 0; M ← MI ; do { for j = 1 to NE VW GY ← CalculaPotencial(G, W, B, Y, P ); (cG ) ← Sorteo(G, VW GY ); VerificaClasifica(cG ); end for AplicarCondicionesContorno(G, W, B, cf in ); VW G ← CalculaPotencial(G, W, B, cf in ); Canal ← CanalTrayectoria(VW G ); if (Canal 6= ∅) (Libre, Y ) ← CompruebaClasifica(Canal); if (Libre == T RU E) T rayectoria ← Trayectoria(Canal); return T rayectoria; else A ← ObtencionCeldasA (VW G ); G ← (G − A); end if ++i; if (M < MM ) ++M ; end if else G ← (G ∪ A); end if } while ((Libre == F ALSE) AN D (i < NI )); return T rayectoria; end; Algoritmo 6.5: Algoritmo panificador de trayectorias, PHM avanzado.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

123

Figura 6.13: Secuencia ejemplo del planificador PHM avanzado. Progreso de la secuencia: de izquierda a derecha y de arriba hacia abajo. Las celdas grises son de contenido desconocido, las negras son de obst´aculo, las blancas son de espacio libre, las azules son las celdas que se excluyen del sorteo de muestreo y las rojas son las celdas de la trayectoria.

124

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

Figura 6.14: Secuencia determinista de muestreo en 2D. El n´ umero indica el orden de la secuencia.

En la figura 6.13 se muestra un ejemplo de ejecuci´on de la implementaci´on del planificador PHM avanzado (la lectura de la secuencia de las figuras se realiza de izquierda a derecha y de arriba a abajo)

6.4.

Exploraci´ on basada en la detecci´ on de colisi´ on: Planificador PHM color

En el planificador de trayectorias PHM se necesita disponer de una funci´on detectora de colisi´on que devuelva la distancia de colisi´on tanto positiva, si el punto de configuraci´on muestreado est´a en el espacio libre, como negativa (penetraci´on) si el punto cae en el espacio obst´aculo. La obtenci´on de una funci´on de este tipo entra˜ na una dificultad no resuelta actualmente, por eso se propone un m´etodo alternativo basado en funciones detectoras de colisi´on que no devuelven la distancia de colisi´on. En este apartado se considera un planificador que realiza una clasificaci´on probabil´ıstica de celdas: Una determinada celda se considera libre o de obst´aculo, con una probabilidad dada, en funci´on de la existencia de colisi´on o no en una serie de puntos obtenidos por muestreo y pertenecientes a dicha celda. Para un n´ umero limitado de puntos de muestreo del subespacio cubierto por una celda, se obtiene una cobertura m´as uniforme de muestreo utilizando alg´ un m´etodo determinista de muestreo, que realizando un muestreo aleatorio. Por este motivo, en el subapartado siguiente, se propone un procedimiento que sirve para generar, secuencialmente, puntos pertenecientes a cualquier celda; de forma que cada punto estar´a separado lo m´as posible de sus predecesores y ser´a coincidente con el punto central de una subcelda descendiente. Este m´etodo es presentado por I˜ niguez y Rosell en la publicaci´on [4] del apartado 7.1.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

6.4.1.

125

Secuencia determinista de muestreo

La secuencia determinista de muestreo n Sm (k) = {s1 , s2 , s3 , . . . , sk , . . .},

(6.26)

aplicada a una celda cnm , provee una ordenaci´on de c´odigos de celdas cuyos centros son muestras que cubren incrementalmente y uniformemente dicha celda. El elemento gen´erico de la secuencia (sk ), vendr´a expresado por: sk = {(m + h), q(m+h)(k) }

(6.27)

Es decir, sk viene definido por el nivel de rejilla (m + h) y las coordenadas de la celdilla sobre dicha rejilla q(m+h)(k) , de acuerdo con el esquema num´erico de identificaci´on de celdas presentado en el apartado 5.2.1. En el ejemplo de la figura 6.14, las primeras a muestrear de la celda c20 , que contiene todo el C-espacio en 2D, son: S02 (k) = {{1, (0, 0)}, {1, (1, 1)}, {1, (0, 1)}, {1, (1, 0)}, {2, (0, 0)}, {2, (2, 2)}, {2, (2, 0)}, }, {2, (0, 2)}, {2, (1, 1)}, {2, (3, 3)}, . . . } (6.28) Y la secuencia generada para muestrear la celda etiquetada con {1, (1, 0)} es: S12 (k) = {{2, (2, 0)}, {2, (3, 1)}, {2, (3, 0)}, {2, (2, 1)}, . . . }

6.4.2.

(6.29)

Clasificaci´ on probabil´ıstica de celdas

La clasificaci´on de una m-celda se realiza dependiendo del resultado de la comprobaci´on de colisiones del punto configuraci´on de su centro y de los centros de sus subceldas. Las configuraciones son consecutivamente obtenidas n por la secuencia determinista de muestreo Sm (k), donde el m´aximo n´ umero de configuraciones generada para la clasificaci´on de celdas es: J=

M X

2n(i−m)

(6.30)

i=m

que se corresponde con el n´ umero de celdas compuesto por la m-celda y todas sus subceldas de niveles comprendidos entre m + 1 y M . La clasificaci´on de una celda se realiza habiendo generado y evaluado solamente j < J configuraciones, y esta clasificaci´on es actualizada cada vez que se toman nuevas muestras. Suponiendo que se han generado j configuraciones y que se han detectado libres, entonces la celda puede clasificarse como libre con un cierto grado de certidumbre que se denominar´a color de certidumbre y que se define como: Ccolor =

j J

(6.31)

Puede realizarse una clasificaci´on an´aloga si todas las j configuraciones son detectadas como obst´aculo. En el ejemplo de la figura 6.15a se muestra la

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

126

y

y

a

15

3

14

11 3

2 6

7

2

10

W

0 13

1

9

16

1 5

0

12 Ccolor=17/21

4 8

0

1

3

2

x

y

x y

b 3

2 6

7

B

0 1

Ccolor=8/21

4

5 x y

x

y

c 3

2

7

G

W

G

Ccolor=2/5

B

B

Ccolor=2/5

Ccolor=1/5

B

6 0 1

G

4

5 x

x

Figura 6.15: Clasificaci´on de celdas: a) Celda libre. b) Celda obst´aculo. c) Celda subdividida en libre, obst´aculo y subceldas desconocidas.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

127

clasificaci´on como celda libre W con (Ccolor = 17/21) y la figura 6.15b muestra la clasificaci´on de una celda obst´aculo B con (Ccolor = 8/21). Cuando todas las J configuraciones son comprobadas como libres de colisi´on u obst´aculo, el color de la celda es conocido con certeza (hasta el m´aximo nivel de resoluci´on), entonces (Ccolor = 1). Las celdas que todav´ıa no han sido exploradas son agrupadas en G y les corresponde (Ccolor = 0). Cuando algunas de las j configuraciones son comprobadas como libres y otras como obst´aculo, la celda no puede ser clasificada como libre ni como obst´aculo puesto que contiene configuraciones libres y configuraciones obst´aculo. En este caso la celda se subdivide en subceldas y as´ı consecutivamente hasta obtener celdas con configuraciones solamente libres o solamente obst´aculo, o sin ninguna muestra. En el primero de los dos casos las subceldas son clasificadas como libres u obst´aculo (con su propio grado de certidumbre), y las subceldas vac´ıas son clasificadas como subceldas desconocidas. Como ejemplo, la figura 6.15c muestra una celda donde las siete primeras muestras se clasifican como obst´aculo pero la octava (muestra n´ umero 7) se clasifica como libre. La celda se particiona en: una M -celda blanca con (Ccolor = 1), dos celdas negras con (Ccolor = 2/5), una celda negra con (Ccolor = 1/5) y tres M -celdas grises.

6.4.3.

Muestreo de celdas

El muestreo de celdas se realiza con el objetivo de explorar el C-espacio desconocido y encontrar una canal a trav´es de las celdas libres, desde una celda bini a una celda destino bf in . Esta exploraci´on de celdas debe de hacerse de forma que sea posible encontrar soluci´on habiendo explorado un n´ umero de celdas lo menor posible. Con este objetivo se establece una probabilidad de muestreo de una celda ck condicionada por el peso ω(ck ): ω(ck ) = ω1 (ck ) · ω2 (ck ) · ω3 (ck )

(6.32)

Estos pesos tienen en cuenta los tres items siguientes: La probabilidad de muestrear una celda incrementa con su tama˜ no, tal como se propuso en la secci´on 6.11, ecuaci´on (6.4). Es decir, si ck es una m-celda entonces: ω1 (ck ) = 2−mn (6.33) Utilizando la ecuaci´on (6.24): ω2 (ck ) =

Vcontorno − VW GY (ck ) , Vcontorno − Vapoyo

(6.34)

La probabilidad de muestrear una celda debe incrementarse con la incertidumbre de color: ω3 (ck ) = 1 − Ccolor

(6.35)

128

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

Figura 6.16: Ejemplo del planificador PHM color, secuencia de izquierda a derecha y de arriba abajo: Las celdas verdes no tienen ning´ un punto de muestreo. Las blancas son aquellas en las que se ha comprobado todos sus posibles puntos de muestreo y han resultado libres, mientras que ser´ıan negras si hubiesen resultado de obst´aculo. Las celdas gris claro son aquellas en las que se han comprobado un n´ umero de puntos m´as peque˜ no de los posibles y han resultado libres, mientras que si hubiesen resultado de obst´aculo ser´ıan gris oscuro. Las celdas con puntos libres y de obst´aculo se subdividen. Y las celdas rojas corresponden al canal de trayectoria.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

129

Las celdas grises tienen ω3 (ck ) = 1 y las celdas cuya certidumbre de color es Ccolor = 1 tienen ω3 (ck ) = 0 y no ser´an seleccionadas para el muestreo. La funci´on integral para el sorteo de celdas, se obtiene como en la ecuaci´on (6.5): F (ck ) =

k X i=1

ω(ci ) ∀ci ∈ G y

1 ≤ k ≤ NG

(6.36)

Finalmente, el muestreo se realiza utilizando la funci´on integral F (ck ) definida sobre las celdas que entran en el sorteo y siguiendo el procedimiento explicado en el apartado (6.1.1).

6.4.4.

Algoritmo del planificador PHM color

El algoritmo de planificaci´on explora iterativamente el C-espacio e intenta encontrar una trayectoria usando las funciones arm´onicas como guia. Para ello dispone de las etapas siguientes: Muestrea un n´ umero de celdas del C-espacio siguiendo el procedimiento del apartado 6.4.3 y la selecci´on de un punto de muestreo, dentro de la celda, siguiendo la secuencia determinista del apartado 6.4.1. Clasifica las celdas nuestreadas siguiendo el procedimiento detallado en el apartado 6.4.2. De acuerdo con los resultados obtenidos en el muestreo de una celda, se actualiza la certidumbre de color o se particiona en subceldas. Computa la funci´on arm´onica sobre las celdas grises y blancas, fijando las celdas negras a un valor alto y la celda destino a un valor bajo. Encuentra un canal de celdas entre cini y cf in siguiendo el gradiente o m´axima pendiente. Si el canal no existe vuelve al paso 1. El canal soluci´on encontrado se compone de celdas, en general de diferentes tama˜ nos, blancas (cada una con su color de incertidumbre) y/o celdas grises (totalmente inexploradas). Entonces: 1. Mientras existan celdas en el canal de tama˜ no superior a la m´axima resoluci´on (M -celdas), se cambia el dominio del C-espacio por el subdominio definido por el canal encontrado. 2. Se llama recursivamente a la funci´on principal del algoritmo, volviendo al paso anterior. En el ejemplo de la figura 6.16 se muestra un C-espacio 2D en diferentes fases del algoritmo planificador de trayectorias leidas de izquierda a derecha y de arriba a abajo. La primera figura muestra la situaci´on inicial con una discretizaci´on no regular cubriendo todo el C-espacio, efectuada sobre la disposici´on conocida de cini y cf in . La u ´ltima figura muestra el canal soluci´on

130

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

obtenido que enlaza la celda inicial y final. PlanificadorTrayectoriaColor(cini ,cf in ) begin (cini , cf in ) ← MarcaInicioFin; G ← DiscreInicial (cini , cf in ); B ← ∅; W ← ∅; Y ← {cini , cf in }; T rayectoria ← ∅; i ← 0; M ← MI ; do { for j = 1 to NE VW GY ← CalculaPotencial(G, W, B, Y ); (cG ) ← Sorteo(G, VW GY ); ChequeoDeterministaClasifica(cG ); end for AplicarCondicionesContorno(G, W, B, cf in ); VW G ← CalculaPotencial(G, W, B, cf in ); Canal ← CanalTrayectoria(VW G ); if (Canal 6= ∅) (Libre) ← CompruebaResolucionCanal(Canal); if (Libre == T RU E) T rayectoria ← Trayectoria(Canal); return T rayectoria; else if PlanificadorTrayectoriaColor(cini ,cf in ); end if end if } while ((Libre == F ALSE) AN D (i < NI )); return T rayectoria; end;

Algoritmo 6.6: Algoritmo panificador de trayectorias PHM color.

En el algoritmo 6.6 se muestra el procedimiento detallado del planificador PHM color, donde las nuevas funciones utilizadas son: ChequeoDeterministaClasifica(cG ): siguiendo la secuencia determinista, presentada en el apartado 6.4.1, comprueba un punto dentro de una celda. De acuerodo con el resultado actualiza la variable color y, si procede, subdivide la celda. CompruebaResolucionCanal(Canal): comprueba la resoluci´on de todas las celdas que componen el canal, si son todas de la m´axima resoluci´on devuelve TRUE y sino devuelve FALSE.

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

6.5.

131

Reducci´ on del tiempo de computaci´ on respecto de una discretizaci´ on uniforme

La base del m´etodo propueto en esta tesis es la descomposici´on aleatoria de celdas en 2n -tree; obteni´endose una discretizaci´on parcial del espacio de configuraciones sobre la que se calcula, en multirresoluci´on jer´arquica, una funci´on potencial arm´onica, subarm´onica o superarm´onica. Tal como se aprecia en la figura 6.1, la reducci´on de celdas de discretizaci´on (que ser´ıan puntos de c´alculo) respecto de una discretizaci´on regular, es tanto mayor cuanto menor es el contorno total de obst´aculo, es decir, el espacio libre y el espacio de obst´aculo est´an agrupados (figura 5.9) Por otro lado, cuanto mayor es el n´ umero de grados de libertad, mayor ser´a el n´ umero de celdillas de m´axima resoluci´on, contenidas en una celda de un determinado nivel de partici´on m (en un espacio de n dimensiones y nivel m´aximo de partici´on M , una celda de nivel de partici´on m contine 2(M −m)n celdillas de m´axima resoluci´on). Por consiguiente, la reducci´on de puntos de c´alculo en el espacio libre (al no particionar celdas que no contengan obst´aculo) y en el espacio obst´aculo (al no particionar celdas contenidas completamente en un obst´aculo) ser´a mayor cuanto mayor sea el n´ umero de grados de libertad. Esto es as´ı porque el contorno es el u ´nico que requiere la m´axima discretizaci´on, y este siempre es de una dimensi´on menos que el espacio al que pertenece. En contrapartida a la importante reducci´on de puntos de c´alculo, se requiere un m´etodo de identificaci´on de celdas que debe ser lo m´as sencillo posible para que el precio pagado no sea excesivo. Este requisito lo cumple el m´etodo propuesto en esta tesis en el subapartado 5.2.1.

6.6.

Aportaci´ on

En este cap´ıtulo se presenta el planificador de trayectorias basado en la combinaci´on del uso de las funciones arm´onicas, calculadas sobre una multirresoluci´on jer´arquica adaptativa, y la aplicaci´on de m´etodos aleatorios guiados. Varias son las aportaciones que, para la obtenci´on de este planificador, se proponen en este cap´ıtulo: - Realizaci´on del estudio probabil´ıstico del m´etodo de muestreo propuesto utilizando un modelo poblacional de celdas. - Obtenci´on de un m´etodo probabil´ıstico de muestreo de celdas, dirijido y condicionado por una funci´on potencial arm´onica. En este caso se presupone que la funci´on detectora de colisi´on utilizada reporta la distancia de colisi´on. - Propuesta de otro planificador aplicable para el caso com´ un de utilizar funciones detectoras de colisi´on que no reportan distancia de colisi´on. En este m´etodo se plantea la exploraci´on probabil´ıstica de celdas, tambi´en guiada por una funci´on potencial, pero en este caso cada vez que se selecciona una

132

´ CAP. 6. FUNCIONES ARMONICAS Y MUESTREO

celda se mustrea un punto de la misma, siguiendo una secuencia determinista. Cada celda se clasifica y procesa dependiendo de la relaci´on de puntos, libres y de obst´aculos obtenidos.

Cap´ıtulo 7 Conclusiones Las distintas aportaciones de cada cap´ıtulo a la tesis han sido referenciadas al final de los mismos; as´ı que en este cap´ıtulo, dedicado a las conclusiones, se destacan las consideradas m´as relevantes. En la figura 7.1 se reproduce el organigrama del planteamiento de objetivos y desarrollo de la tesis, as´ı como la situaci´on de los cap´ıtulos en dicho contexto, para centrar las aportaciones principales obtenidas: 1. Como primera propuesta fundamental, en el cap´ıtulo 3, se obtienen alternativas compensatorias de la distribuci´on del gradiente, basadas en la utilizaci´on de funciones subarm´onicas y superarm´onicas, que reduce los inconvenientes de los m´etodos cl´asicos de planificaci´on de movimientos basados en funciones arm´onicas. 2. En el cap´ıtulo 4 se analiza la funci´on potencial, dependiente del tiempo, como soluci´on de un modelo que incluye un obst´aculo m´ovil, en base al cual se dise˜ na una estrategia de control correctivo de posici´on y de movimientos en entornos con obst´aculos m´oviles. 3. En el cap´ıtulo 5, como tercera propuesta fundamental, se plantea una discretizaci´on jer´arquica adaptativa en multirresoluci´on y por fases, que reduce dr´asticamente el n´ umero de puntos de c´alculo requeridos en la soluci´on num´erica de las funciones arm´onicas, posibilitando su uso para problemas con un n´ umero de grados de libertad mayor que los resueltos habitualmente mediante estos m´etodos. 4. Como cuarta y u ´ltima propuesta fundamental, en el cap´ıtulo 6, se propone la combinaci´on de t´ecnicas de muestreo con el c´alculo de funciones arm´onicas mediante un m´etodo de exploraci´on aleatorio guiado (PHM) aplicado a un espacio de configuraciones discretizado jer´arquicamente sobre el que se va recalculando la funci´on arm´onica, permitiendo encontrar soluciones de forma m´as r´apida al centrarse solamente en la parte del espacio relevante al problema. 133

134

CAP. 7. CONCLUSIONES

Figura 7.1: Organigrama del desarrollo de la tesis.

7.1.

Publicaciones realizadas con el desarrollo de la tesis

Las publicaciones afines, m´as representativas, realizadas durante el desarrollo de la tesis son las siguientes: [1] Rosell J. and I˜ niguez P. (2002) A Hieralchical and Dynamic Method to Compute Harmonic Functions for Constrained Motion Planning. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS’02, pp. 2334-2340. En este trabajo se introduce el c´alculo de las funciones arm´onicas en un espacio de configuraciones discretizado jer´arquicamente, siguiendo el m´etodo presentado en el apartado 5.2. [2] I˜ niguez P. and Rosell J. (2003) Probabilistic Harmonic-function-based Method for Robot Motion Planning. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS’03. En esta publicaci´on se propone el planificador Probabilistic Harmonic-functionbased Method (PHM), el cual plantea la intercalaci´on de m´etodos de muestreo aleatorio y de c´alculo de funciones arm´onicas, tal como se describe en el apartado 6.1. [3] I˜ niguez P. and Rosell J. (2003) Combining Harmonic Functions and Random Sampling in Robot Motion Planning: A performance Evaluation. Proceedings of the IEEE Symposium on Assembly and Task Planning, pp. 253-258. Aqu´ı se realiza una evaluaci´on te´orica del m´etodo anterior, obteniendo la relaci´on entre la probabilidad de fallo y el n´ umero de muestras, tal como se

CAP. 7. CONCLUSIONES

135

detalla en el apartado 6.1.4. [4] P. I˜ niguez and J. Rosell (2005) Combining Harmonic Functions and Random Sampling in Robot Motion Planning: A lazy approach. Proceedings of the IEEE Symposium on Assembly and Task Planning - ISATP 2005. Fundamentalmente, este trabajo mejora al anterior realizando un control pasivo y algor´ıtmico de la resoluci´on m´axima aplicable a la discretizaci´on jer´arquica, mientras intenta encontrar una trayectoria, y controlando las regiones donde focalizar el inter´es, mejorando as´ı la eficiencia computacional, estas mejoras est´an incluidas en los apartados 6.2 y 6.3. [5] J. Rosell and P. I˜ niguez (2005) Path Planning using Harmonic Functions and Probabilistic Cell Decomposition. Proceedings of the IEEE Int. Conf. on Robotics and Automation - ICRA 2005, pp. 1815-1820. Best Paper Award Finalist. Este trabajo presenta la extensi´on del planificador PHM considerando el caso de detectores de colisi´on sin informaci´on de distancia. Mediante el uso de muestreo determinista se eval´ ua la probabilidad de que una celda sea libre o de colisi´on, y se adapta el m´etodo a esta incertidumbre. Este enfoque es desarrollado en el apartado 6.4. [6] I˜ niguez P. and Rosell J. (2009) Path planning using sub- and superharmonic functions. Proceedings of the 40th International Symposium on Robotics, ISR’09 , pp.319-324. En este trabajo se plantea la utilizaci´on de las funciones superarm´onicas para simular los obst´aculo y las subarm´onicas para generar un pozo de potencial en un entorno pr´oximo al punto destino. De esta forma la funciones potenciales obtenidas presentar´an muchas menos zonas planas que las obtenidas con las funciones arm´onicas. En el apartado 3.5 se analiza esta propuesta.

7.2.

Trabajos futuros

Dentro de esta tesis, varias son las l´ıneas sobre las que se podr´ıa incidir con la finalidad de desarrollar trabajos que permitan mejorar el planificador de trayectorias propuesto: Incorporaci´on de las propuestas en entornos de planificaci´on que incluyan el modelado geom´etrico de problemas reales y posible adaptaci´on para el caso de disponer solo de distancias de colisi´on positivas (no de penetraci´on). Aplicaci´on a problemas de rob´otica m´ovil mediante la adaptaci´on del m´etodo a situaciones en las que el entorno solo se conoce parcialmente y del que se puede recabar m´as informaci´on mediante los sensores del robot.

136

CAP. 7. CONCLUSIONES

Ap´ endice A Demostraciones Con la finalidad de obtener una mayor claridad de exposici´on, en este ap´endice se sit´ uan las demostraciones matem´aticas de las proposiciones, teoremas y ecuaciones m´as relevantes, utilizadas en el desarrollo de la tesis. Sin embargo, a pesar de presentarse en el ap´edice, el contenido de los apartados A.1, A.2 y A.3 son tambi´en aportaciones b´asicas obtenidas en la consecuci´on de los ojetivos de la tesis.

A.1.

Principio de localidad de una funci´ on

Proposici´ on para dos dimensiones: Dada una funci´on f (x, y), definida en un dominio Ω ⊂ R2 y de clase C 2 . Dado un disco de radio R centrado en (x0 , y0 ): D = {(x, y) ∈ Ω|kx − x0 , y − y0 k < R}. Entonces se cumple que:

1 f (x0 , y0 ) = 2πR

I

L

f (x, y)dl −

R2 div (grad f )0 4

(A.1)

donde la integral de l´ınea se realiza sobre la circunferencia contorno del disco. Prueba: Aplicando Taylor en el punto P0 de la superficie z = f (x, y), de clase C 2 : 1 ∂f 1 ∂f (x − x0 ) + (y − y0 )+ f (x, y) = f (x0 , y0 ) + 1! ∂x (x0 ,y0 ) 1! ∂y (x0 ,y0 ) 1 ∂ 2 f 1 ∂ 2 f 2 (x − x0 ) + (y − y0 )2 + 2! ∂x2 (x0 ,y0 ) 2! ∂y 2 (x0 ,y0 ) 2 ∂ 2 f (x − x0 )(y − y0 ) (A.2) 2! ∂x∂y (x0 ,y0 ) 137

´ APENDICE

138

Figura A.1: Principio de localidad de una funci´on. Integrando en torno de la circunferencia (x0 +Rcosθ, y0 +Rsenθ, 0), teniendo en cuenta que dl = Rdθ (figura A.1): Z Z I Z 2π 1 ∂f 2π 1 ∂f 2π f (x, y)dl = f (x0 , y0 ) RcosθRdθ+ RsenθRdθ+ Rdθ+ 1! ∂x 0 1! ∂y 0 L 0 Z Z 1 ∂ 2 f 2π 2 2 1 ∂ 2 f 2π 2 R cos θRdθ + R sen2 θRdθ+ 2! ∂x2 0 2! ∂y 2 0 Z 2π 2 ∂ 2f R2 cosθsenθRdθ (A.3) 2! ∂x∂y 0 Los t´erminos correspondientes a la derivadas impares y cruzadas se anulan:   I 1 ∂ 2f ∂ 2f f (x, y)dl = f (x0 , y0 )2πR + + 2 πR3 , (A.4) 2 2! ∂x ∂y L despejando se obtiene la expresi´on A.1. Proposici´ on para n dimensiones: Sea u = f (q1 , q2 , . . . , qn ) una funci´on de n variables, definida en un dominio Ω en Rn e infitamente derivable. Sea una bola B ∈ Rn , centrada en r0 = (q10 , q20 , . . . , qn0 ) y de radio R. El valor de la funci´on en r0 vendr´a dada por la expresi´on 1 f (r0 ) = Sn

I

S

f (r0 + Rv)dS −

R2 2 R4 4 ∇ f (r0 ) − ∇ f (r0 ) − . . . 2 · 2! 8 · 4!

(A.5)

Prueba: Aplicando Taylor en el punto r0 = (q10 , q20 , . . . , qn0 ) siguiendo la direcci´on del vector unitario v (figura A.2) f (r0 + tv)

=

f (r0 ) +

1 ∂ 2f 2 1 ∂ 3f 3 1 ∂f t + t + t + . . . (A.6) 1! ∂v 2! ∂v2 3! ∂v3

´ APENDICE

139

Figura A.2: Hiperesfera centrada en el punto P0

donde t es un escalar y las derivadas direccionales son sobre la direcci´on se˜ nalada por el vector v = (h1 , h2 , . . . , hn ). Integrando en una hiperesfera centrada en el punto, de radio t = R I

R f (r0 + Rv)dS = f (r0 ) dS + 1! S S I

I

S

I 2 ∂f R2 ∂ f dS + dS+ ∂v 2! S ∂v2 I 3 ∂ f R3 dS + . . . (A.7) 3! S ∂v3

Las derivadas direccionales deben permanecer dentro de la integral, puesto que el vector v apunta al diferencial dS sobre la superficie de integraci´on (figura A.2) Utilizando la expresi´on de la derivada direccional: n

X ∂f ∂f = ∇u · v = hi ∂v ∂q i i=1

(A.8)

y aplicando la regla de la cadena n X ∂ 2f ∂ 2f = hi hj ∂v2 i,j=1 ∂qi ∂qj n X ∂ 3f ∂ 3f = hi hj hk ∂v3 i,j,k=1 ∂qi ∂qj ∂qk

(A.9)

(A.10)

´ APENDICE

140

se obtiene: I

n

I n R2 X ∂ 2 f hi dS+ hi hj dS+ 2! i,j=1 ∂qi ∂qj S S I n ∂ 3f R3 X hi hj hk dS + . . . (A.11) 3! i,j,k=1 ∂qi ∂qj ∂qk S

R X ∂f f (r0 +Rv)dS = f (r0 ) dS+ 1! i=1 ∂qi S S I

I

Los t´erminos correspondientes a la derivadas impares y cruzadas se anulan, mientras que permanecen los correspondientes a las derivadas pares: Aplicando el teorema del valor medio a la parte integral del primer sumatorio Z I I S π 1 · cos θdθ = 0, (A.12) hi dS = hiM ED dS = hiM ED S = π 0 S S y en la parte integral del segundo sumatorio, si i 6= j, se cumple I hi hj dS = S(hi hj )M ED = ShiM ED hjM ED = 0

(A.13)

S

ya que qi y qj son ortogonales, y si i = j I I Z S π S 2 2 hi hj dS = hi dS = ShiM ED = cos2 θdθ = π 0 2 S S y I

S

h4i dS

=

Sh4iM ED

S = π

Z

π

cos4 θdθ = 0

S 8

(A.14)

(A.15)

Similarmente ocurre con los dem´as t´erminos, con lo cual la expresi´on (A.11) se simplifica I

I

n

n

R4 X ∂ 4 f R2 X ∂ 2 f + S + ... f (r0 + Rv)dS = f (r0 )S + S 2 · 2! i=1 ∂qi2 8 · 4! i=1 ∂qi4 S

(A.16)

f (r0 + Rv)dS = f (r0 )S + S S

R2 2 R4 4 ∇u + S ∇ u + . . . (A.17) 2 · 2! 8 · 4!

y, aislando f (r0 ) se obtiene (A.5).

A.2.

´ Area e integraci´ on de una hiperesfera

La ecuaci´on de una hiperesfera en n dimensiones, centrada en el origen, es: x21 + x22 + x23 + · · · + x2n = r2

(A.18)

´ APENDICE

141

y, el vector unitario normal: v=

x2 x3 xn x1 x 1 + x 2 + x3 + · · · + x n r r r r

Cuya divergencia es:  div v = ∂x∂ 1 x1 + ∂x∂ 2 x2 +

∂ x ∂x3 3

+ ··· +

∂ x ∂xn n

div v =

 ·

x1 x r 1

n r

se obtiene:

n r

R

V

dV = 1 ·

H

S

x2 x r 2

+

x3 x r 3

+ ··· +

(A.20)

Aplicando el teorema de la divergencia: Z I (∇ · v) dV = v · dS V

+

(A.19)

(A.21)

S

⇒

dS

n V n = Sn r

(A.22)

que relaciona el volumen (Vn ) con la superficie (Sn ) de una hiperesfera de n dimensiones y radio r. Para el espacio tridimensional, la ecuaci´on de la superficie esf´erica es: x2 +y 2 +z 2 = r2 , y tomando x como par´ametro, se obtiene: y 2 +z 2 = r2 −x2 , que es la superficie de un c´ırculo. Por lo tanto el volumen del cilindro difrencial es: dV3 = π(r2 − x2 )dx ,

(A.23)

y el volumen de la esfera, ser´a: V3 = π

Z

+r −r

(r2 − x2 )dx ,

donde realizando el cambio: x = r cos ϑ, se obtiene: Z π 3 sin3 ϑdϑ V3 = πr

(A.24)

(A.25)

0

Esta integral se resuelve utilizando el cambio: t = cos ϑ, obt´eniendose: 4 V3 = πr3 , 3

(A.26)

xn x r n




´ APENDICE

142 y

3 S3 = V3 = 4πr2 r

(A.27)

Siguiendo el mismo procedimiento para n = 4, se tiene: 1o - Ecuaci´on de la hiperesfera: x21 + x22 + x23 + x24 = r2 2o - Parametrizando respecto a x1 , tenemos la esfera: x22 + x23 + x24 = r2 − x21 3o - Hipercapa esf´erica diferencial: dV4 = 4o - Volumen de la hiperesfera 4d: V4 =

4π 2 (r 3

4π 3

R +r −r

3

− x21 ) 2 dx1 3

(r2 − x21 ) 2 dx1

5o - Realizando el mismo cambio que antes (x1 = r cos ϑ), se obtiene: Z 4 4 π 4 V4 = πr sin ϑdϑ (A.28) 3 0 6o - Volumen de la hiperesfera: V4 =

π2 4 r 2

(A.29)

7o - Superficie de la hiperesfera: 4 S4 = V4 = 2π 2 r3 r Y repitiendo el procedimiento para n = 5, se tiene: Z π 2 r5 π 5 sin ϑdϑ V5 = 2 0 V5 =

8π 2 5 r 15

5 8π 2 4 S5 = V 4 = r r 3

(A.30)

(A.31)

(A.32)

(A.33)

Siguiendo el m´etodo de inducci´on sobre n, se obtienen las expresiones gen´ericas: Z π n sinn ϑdϑ (A.34) Vn = kn r 0

Vn = kn+1 rn

(A.35)

´ APENDICE

143

Figura A.3: Estudio de la convexidad.

Sn = nkn+1 rn−1

(A.36)

Donde kn+1 es: n

kn+1 y la funci´on Γ (x) es: Γ (x) =

Z

π2
= n Γ 2 +1

(A.37)

+∞

tx−1 e−t dt; 0+

x ∈ (0, +∞)

(A.38)

Tambi´en, teniendo en cuenta las ecuaciones (A.34) y (A.22) se tiene: Z π n−1 sinn ϑdϑ (A.39) Sn = nkn r 0

Las expresiones (A.34) y (A.39) parametrizan la integraci´on de una hiperesfera en funci´on de una u ´nica variable.

A.3.

Propiedades de las funciones arm´ onicas, superarm´ onicas y subarm´ onicas

A.3.1.

Concavidad, convexidad de las funciones arm´ onicas y superarm´ onicas

Considerando tres puntos alineados (r0 , r1 y r2 ) en Rn (figura A.3), de forma que r0 tiene una posici´on variable intermedia, cuya expresi´on param´etrica

´ APENDICE

144 se puede definir de dos maneras: r0 = r1 + (r2 − r1 )ϑ r0 = r2 − (r2 − r1 )(1 − ϑ) sumando y simplificando, se tiene

 

, 0