Optimización de misiones interplanetarias mediante el uso de algoritmos genéticos
Abstract
En el diseño de misiones espaciales, ya sea dentro de la influencia terrestre o enfocadas al espacio profundo, es primordial realizar una optimización en pos de minimizar el combustible necesario para lograr el objetivo, haciendo asumible y atractiva económicamente la exploración y explotación del espacio. Dada la no linealidad del problema y la cantidad de variables en juego, es necesario proceder por etapas, y encontrar primero una buena solución bajo ciertas hipótesis para analizarla posteriormente en detalle. Para esta primera etapa, en la que se centra este proyecto, existen distintas alternativas de optimización. La que se va a abordar aquí es el uso de algoritmos genéticos, de aparición en 1970 y fuerte crecimiento en el campo de estudio en la última década gracias a los avances computacionales, ofreciendo grandes ventajas frente a otros métodos. Una herramienta con enorme potencial que se desarrollará y pondrá a prueba en esta memoria.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Optimización de misiones interplanetarias mediante el uso de algoritmos genéticos Autor: Jesús F. Ramírez Sánchez Tutores: Rafael Vázquez Valenzuela Julio Cesar Sánchez Merino Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020
Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Optimización de misiones interplanetarias mediante el uso de algoritmos genéticos Autor: Jesús F. Ramírez Sánchez Tutores: Rafael Vázquez Valenzuela Profesor Titular Julio Cesar Sánchez Merino Personal Investigador en Formación Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020
Trabajo Fin de Grado: Optimización de misiones interplanetarias mediante el uso de algoritmos genéticos Autor: Jesús F. Ramírez Sánchez Tutores: Rafael Vázquez Valenzuela Julio Cesar Sánchez Merino El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Resumen En el diseño de misiones espaciales, ya sea dentro de la influencia terrestre o enfocadas al espacio profundo, es primordial realizar una optimización en pos de minimizar el combustible necesario para lograr el objetivo, haciendo asumible y atractiva económicamente la exploración y explotación del espacio. Dada la no linealidad del problema y la cantidad de variables en juego, es necesario proceder por etapas, y encontrar primero una buena solución bajo ciertas hipótesis para analizarla posteriormente en detalle. Para esta primera etapa, en la que se centra este proyecto, existen distintas alternativas de optimización. La que se va a abordar aquí es el uso de algoritmos genéticos, de aparición en 1970 y fuerte crecimiento en el campo de estudio en la última década gracias a los avances computacionales, ofreciendo grandes ventajas frente a otros métodos. Una herramienta con enorme potencial que se desarrollará y pondrá a prueba en esta memoria. I
Abstract In the design of space missions, wheter they are under the Earth’s influence or focused on deep space, it is essential to solve an optimization problem in order to minimize the fuel consumption needed to reach the objetive thus making affordable and financially attractive the exploration and explotation of the space. Due to the non-linearity of the problem and the amount of variables involved, it is necessary to proceed step by step, to find a first good solution under certain hyphotesis, and then analyzing it in detail. For this first step, the context in which this work is developed, there are different optimization alternatives. The one which is going to be treated here is the use of genetic algorithms, which appeared in 1970 and has strongly grown in the space field in the last decade thanks to the computation developements, offering big advantages compared with other methods. Genetic algorithms are a tool with an enormous potential which is going to be developed and tested in this document. III
XNotación λTLongitud Verdadera. ϖLongitud del periapsis. uArgumento de la latitud. àPrimer Punto de Aries. ÀSol. ⊕Tierra. cVector de elementos orbitales. c0Vector de elementos orbitales en J-2000. ˙ cDerivada de elementos orbitales. MEP Maniobra de Espacio Profundo. UB Límite superior de una variable. LB Límite inferior de una variable. AMatriz de restricciones lineales. bVector de restricciones lineales. X Número de individuos que conforman la población base. ENúmero de individuos catalogados de élite. NNúmero máximo de generaciones permitidas. tNúmero de integrantes de torneos. ncNúmero de cruces aleatorios por generación. ηc Parámetro de control de similitud descendenciaprogenitores. σc Divisor en la varianza de la distribución Normal en mutación. pmut Probabilidad de un individuo de mutar. pmProbabilidad de mutar de gen en un individuo mutante. pc Probabilidad de un individuo de convertirse en progenitor. TOL Tolerancia de discrepancia entre el valor medio en la población en la Función Objetivo y el del mejor individuo. ∆VTImpulso total empleado en la misión. TVTiempo de Vuelo Total de la misión. JD Días Julianos.
1 Introducción La exploración del espacio profundo ha sido desde los orígenes de la carrera espacial, en el año 1957, un objetivo para muchos países y compañías, pero al alcance de pocos dados sus altos costes y tecnologías requeridas. Sin embargo, en las últimas décadas se puede observar que cada día son más las misiones que se lanzan hacia la exploración de nuestro Sistema Solar, e incluso más allá. En las Figuras 1.1 y 1.2, pueden apreciarse algunos ejemplos en lo relativo a la exploración del Sistema Solar, el autor de estas imágenes expone que "These days, it’s hard to keep track, because there are so many." ("En estos días, es difícil seguir la pista, porque son demasiadas") [ 7 ]. Sin más que mirar los colores de los gráficos puede comprobarse el aumento del número de este tipo de misiones. Figura 1.1 Misiones de espacio profundo de mayor relevancia dentro del Sistema Solar, anteriores al año 2003. Los motivos de este incremento son varios: •Nuevos intereses en la exploración del espacio profundo tales como: –Misiones científicas y de exploración. –Aplicaciones tecnológicas y comerciales. –Misión de colonización de Marte. 1
2Capítulo 1. Introducción • Avances en la tecnología y conocimientos que permiten llegar más lejos y reduciendo costes. • Avances en el cálculo computacional, que permiten la creación de herramientas más potentes para la optimización de trayectorias, con la consecuente reducción de precios y tiempos. Figura 1.2 Misiones de espacio profundo de mayor relevancia dentro del Sistema Solar, posteriores al año 2003 y proyectadas hasta 2018. La optimización de las misiones supone un ahorro considerable para las agencias y la posibilidad de aumentar la carga de pago gracias a la disminución de masa en combustible, con el consiguiente mayor retorno científico. Un proceso obligatorio para cualquier misión, un paso que puede realizarse sin costes materiales pero que conlleva una gran complejidad que exige dividir este procedimiento en varios subprocesos, desde las soluciones más simples hasta su último análisis en detalle y definición de la misión final. El primer subproceso consiste en determinar en función de una serie de objetivos un primer punto de partida, un óptimo global encontrado utilizando hipótesis simplificativas que en posteriores análisis de detalle se eliminan. Este subproceso no es trivial, y envuelve una serie de complejidades y decisiones a tomar que son las abordadas en este trabajo. 1.1 Alcance del proyecto/trabajo Este proyecto tiene como objetivo probar la eficacia del algoritmo de optimización metaheurístico conocido como algoritmo genético en la resolución de un problema de optimización de transferencia interplanetaria con varias Maniobras de Espacio Profundo; e investigar cómo afectan estas maniobras a los trayectos de estas características. El algoritmo en sí mismo, junto a todos los programas que intervienen en el problema, serán creados desde cero en el entorno de MATLAB (a excepción de algunas funciones básicas como @fsolve, y del denominado refinador @fmincon). Para investigar el funcionamiento del algoritmo y la influencia de las MEP (Maniobras de Espacio Profundo) se realizan una serie de simulaciones de misiones a distintos planetas en el plano bidimensional, y a posterior y como ampliación del trabajo una más a Marte con la consideración de problema tridimensional. Para cerrar la investigación, se estudia la influencia de la fecha de salida
1.2 Estructura del documento 3 en la misión y se comparan los resultados con otro autor. Todos los conceptos mencionados quedan convenientemente explicados a lo largo de la memoria. 1.2 Estructura del documento Tras esta breve introducción, el documento se divide en capítulos que abordan cada una de las partes involucradas en el proceso de optimización, los resultados, y las conclusiones obtenidas tras la realización de este trabajo. En primer lugar se tratarán el Modelo e hipótesis (Capítulo 2) utilizados para definir el problema de Mecánica Orbital asociado a la misión a optimizar, incluyendo una breve introducción a esta ciencia, y los conceptos básicos necesarios para entender el problema. En segunda instancia, se expone la herramienta utilizada para llevar a cabo dicha optimización: el Algoritmo genético (Capítulo 3). En este capítulo se incluye de manera similar al anterior una introducción acompañada de los conceptos básicos de la metodología detrás del algoritmo; y se desarrolla dicho algoritmo enfocado al caso que ocupa. En cuarto lugar, se realiza un Análisis de procedimientos del algoritmo genético (Capítulo 4) con el objetivo de ajustar el programa y evaluar el camino seguido por el algoritmo para encontrar el óptimo aplicado a un escenario concreto, determinando así su convergencia y eficacia. Se procede tras el desarrollo del algoritmo a presentar las Simulaciones y resultados (Capítulo 5). Se presentan en esta sección distintos escenarios y los resultados ofrecidos tanto por el algoritmo creado, como los obtenidos mediante otros algoritmos de optimización de uso genérico presentes en MATLAB , para mostrar así una comparativa que valide y demuestre la eficacia de la herramienta desarrollada. Por último, en lo relativo a las partes principales del documento, se presentan las Conclusiones y trabajo futuro (Capítulo 6), que valoran la solución obtenida y realizan propuestas para proseguir el trabajo en un futuro. Adicionalmente, se añaden a la memoria dos apéndices con información que se extrae del documento base para evitar distraer la atención del lector del hilo principal del proyecto. El primero de ellos: Código Algoritmo genético (Apéndice A) reúne el código completo del optimizador desarrollado en MATLAB . El segundo: Herramientas para modelización del problema (Apéndice B) lista todos los subprogramas desarrollados para la resolución del problema, y entra en detalle en aquellos que se consideran más importantes.
2 Modelo e hipótesis El problema que nos atañe consiste no en una simple transferencia entre órbitas en el Sistema Solar, sino en una transferencia orbital entre las posiciones de distintos planetas. A priori pudiera parecer que el problema es similar, pero conlleva una serie de peculiaridades que son tratadas en este capítulo. Comenzando por introducir al lector en los conceptos genéricos básicos en los que se apoya la Mecánica Orbital, los necesarios para entender el problema, se continúa después por recorrer los distintos matices que caracterizan a este, finalizando con las hipótesis, algunas necesarias y otras simplificativas, que se aplican en el desarrollo del trabajo. 2.1 Conceptos fundamentales de la Mecánica Orbital La rama de la Mecánica Orbital de la que se hace uso en este trabajo es la relativa a la dinámica aplicada a un punto y al uso de órbitas Keplerianas. No se tienen por tanto en cuenta inercias, irregularidades en el campo gravitatorio de los planetas (es decir se consideran estos esferas perfectas), campos gravitatorios de otros cuerpos, u otra serie de dificultades que aparecen a medida que se liberan grados de libertad. Las bases del problema se sientan pues en el uso adecuado de las Leyes de Newton , de la fuerza gravitatoria , y de las Leyes de Kepler . Estas últimas son de especial relevancia dentro de la Mecánica Orbital [12] : 1La órbita de cada planeta es una elipse con el Sol en uno de sus focos. 2 Una línea que una el Sol con un planeta barre áreas iguales en tiempos iguales (Ley de Áreas). 3 El cuadrado del periodo orbital de cada planeta es directamente proporcional al cubo del semieje mayor de la elipse de la órbita del planeta. Estas leyes, si bien enunciadas para el Sistema Solar, son perfectamente aplicables bajo la hipótesis del Problema de los dos cuerpos (explicada en la Sección 2.4) a un satélite, sustituyendo planeta por este, y manteniendo Sol o cambiándolo por el planeta del sistema de dos cuerpos. Son estas leyes deducidas primeramente de la observación del cielo, y obtenibles a partir de la correcta manipulación de las leyes de Newton y la ecuación de la energía; las que permiten determinar el movimiento de un satélite bajo la influencia del campo gravitatorio de un cuerpo principal a través del espacio, las que permiten determinar la cónica en la que se mueve. Es necesario, para el correcto entendimiento de los resultados que se muestran en este documento, explicar cómo se determinan estas cónicas que sigue como trayectoria el satélite, y cómo se relacionan 5
6Capítulo 2. Modelo e hipótesis con su posición y velocidad así como con el tiempo transcurrido en ella. Por ello los siguientes apartados, desde el 2.1.1 hasta el 2.1.3 incluido, indagan en la cuestión, y presentan el conocido como Problema de Lambert , que no es más que el mismo problema con distintas variables de entrada, consistente en encontrar la cónica que une dos posiciones en el espacio para un determinado tiempo de vuelo. 2.1.1 Elementos orbitales, determinación de órbita Obtención de integrales primeras La primera Ley de Kepler puede deducirse para el Problema de dos cuerpos mediante la manipulación matemática de la ecuación del movimiento del problema en cuestión. Esta ecuación del movimiento presenta la forma: ¨ r=−µr r3,(2.1) donde µ=G(m1+m2) es el Parámetro Gravitacional de las Masas Combinadas (siendo G la constante de gravitación universal), r el vector posición que une ambos cuerpos, y ˙ r , ¨ r y r sus derivadas primera con respecto al tiempo, segunda y módulo respectivamente 1 . Aplicando las siguientes operaciones matemáticas en pos de determinar integrales primeras de la ecuación del movimiento, empezando por: [12] ¨ r·˙ r=−µr·˙ r r3,(2.2) r·˙ r=r·v=r˙r−→ ¨ r·˙ r=−µ˙r r2,(2.3) d dt (v2) = d dt (v·v) = 2¨ r·˙ r−→ d dt v2 2=−µ˙r r2=d dt µ r,(2.4) se obtiene la primera de las integrales primeras necesarias para la determinación de la órbita, la ecuación de conservación de la energía específica: v2 2−µ r=cte =ε.(2.5) Procediendo ahora de forma similar desde (2.1), pero haciendo uso de las propiedades del producto vectorial en lugar del escalar esta vez: r∧¨ r=−µr∧r r3=0,(2.6) d dt (r∧˙ r) = ˙ r∧˙ r+r∧¨ r=r∧¨ r−→ d dt (r∧˙ r) = 0.(2.7) Se obtiene con este procedimiento las siguientes tres integrales primeras (es vectorial), la conservación del momento cinético específico h: r∧˙ r=r∧v=h=~ cte.(2.8) Con esta relación es inmediato comprobar la primera parte del enunciado de la primera Ley de Kepler, pues r·h=0 y v·h=0 con h siendo constante, lo que indica que el movimiento se desarrolla en un plano perpendicular a h. 1 Se sigue la notación especificada en el apartado del documento presente en la página VII para todas las variables del desarrollo.
2.1 Conceptos fundamentales de la Mecánica Orbital 7 Para demostrar que este movimiento se corresponde con una trayectoria cónica, es necesaria una última manipulación: ¨ r∧h=−µr∧h r3,(2.9) aplicando la regla del doble producto vectorial a∧(b∧c) = (a·c)b−(a·b)c, se tiene: r∧h=r∧(r∧v) = (r·v)r−(r·r)v=r˙rr−r2v,(2.10) haciendo uso de la igualdad (2.10) en (2.9), y expresando (2.9) como una derivada (donde se ha usado la propiedad de que hes constante) se llega a: d dt (˙ r∧h) = ¨ r∧h+˙ r∧˙ h=¨ r∧h−→ d dt (˙ r∧h) = d dt (v∧h) = −µ˙rr−rv r2=d dt µr r,(2.11) y finalmente, se obtienen las últimas integrales primeras 2 , la conservación del vector excentricidad: v∧h−µr r=~ cte =µe,(2.12) representando e el vector excentricidad ovector de Lagrange. Para terminar, definiendo como θ (la anomalía verdadera) el ángulo formado por r y e , y e (excentricidad) al módulo de e , sin más que multiplicar escalarmente (2.12) por rse llega a: h2−µr=µrecosθ.(2.13) Simplemente reordenando la expresión anterior y despejando r , se obtiene finalmente la ecuación que describe la trayectoria del movimiento, que efectivamente se corresponde con la de una cónica, cerrando este desarrollo de demostración de la primera Ley de Kepler: r= h2 µ 1+ecosθ=p 1+ecosθ.(2.14) Es (2.14) posiblemente una de las ecuaciones más famosas de la Mecánica Orbital, donde p recibe el nombre de parámetro de la cónica. Será precisamente (2.14), sumada a (2.5), (2.8) y (2.12) las que permitan determinar la órbita en la que se mueve el satélite de estudio (o cualquier cuerpo en un problema de dos cuerpos). Cabe mencionar, que la excentricidad ( e ) determina el tipo de cónica, existiendo los siguientes casos: •e=0; Circunferencia. •0<e<1; Elipses. •e=1; Parábola. •e>1; Hipérbolas. Elementos orbitales y sistema de referencia Antes de proceder, llegados a este punto es conveniente introducir los elementos orbitales . Estos son los parámetros mínimos y suficientes para describir sin ambigüedad una órbita en el espacio, y la posición del cuerpo dentro de ella en un instante. Durante este desarrollo ya se han definido tres de ellos: p , e , y θ . Este último lleva implícito el carácter temporal, pues determina en qué punto de la cónica se encuentra el satélite. Dado el origen de estos tres elementos orbitales, son por naturaleza suficientes para definir la cónica de la trayectoria como curva; siendo ahora necesarios para colocarla en el espacio otros tres parámetros. Dicho con otras palabras, los parámetros p y e dibujan la curva en un plano y θ sitúa un punto dentro de esta, quedando por determinar la 2 Puede extrañar el hecho de que aparezcan un total de 7 integrales primeras para un problema de 6 grados de libertad, puede demostrarse que una de estas integrales primeras es combinación lineal del resto: e está contenida en el plano definido por ryv(e·h=0), por lo que solo posee dos grados de libertad.
8Capítulo 2. Modelo e hipótesis orientación de ese plano (conocido como plano/sistema de referencia perifocal ) en el espacio. Pero para poder referir cualquier cosa en un espacio tridimensional, es necesario hablar primero del sistema de referencia. Existen varios sistemas de referencia en la Mecánica Orbital, en función del cuerpo principal que ocupe el problema de estudio (por ejemplo heliocéntrico o geocéntrico) y en función del plano de referencia (perifocal, eclíptico...). Para este trabajo en concreto interesa uno particularmente, que es comentado en la Sección 2.2; pero para introducir los elementos orbitales que faltan se va a utilizar uno genérico que se define en la Figura 2.1: Figura 2.1 Sistema de referencia genérico para órbita. En este sistema de referencia pueden apreciarse el plano de referencia genérico así como el planeta de mayor masa situado en el foco de la cónica (véase el Problema de los dos cuerpos de la Sección 2.4). Se refleja también en la Figura 2.1 el vector de excentricidad e que define la conocida como línea de ápsides, una línea que une el punto de menor distancia entre los dos cuerpos (periapsis) y el foco; el vector h , perpendicular al plano orbital; así como el ángulo θ que siempre se considera positivo en el sentido del movimiento3. Pero en la Figura 2.1 aparecen también nuevos elementos no mencionados con anterioridad, los elementos orbitales restantes. Surge de manera natural un nuevo vector que sirve de herramienta tanto para definir los elementos orbitales, como para pasar de ellos a posición y velocidad, el vector de la línea de nodos ( n ). Este vector tiene como dirección la intersección entre el plano orbital y el de referencia, y sentido hacia el denominado como nodo ascendente, un punto que marca el paso del satélite de la zona inferior del plano de referencia a la superior. 3 Mientras que en elipses se recorre desde 0 º a 360 º , en parábolas, por convención, se hace desde -180 º a 180 º e hipérbolas desde −θ∞aθ∞(obsérvese la Figura 2.4) para evitar ambigüedades en las leyes horarias.
2.1 Conceptos fundamentales de la Mecánica Orbital 9 El primero de los nuevos elementos orbitales es la RAAN ( Ω ,Right Ascension of the Ascending Node), ángulo que forma la línea de nodos con el eje x, comúnmente (y así se tomará en este documento) tal que apunta al primer punto de Aries à . El segundo es el ángulo formado entre el vector excentricidad y n cuya función es situar el periapsis de la órbita, conocido como Argumento de periapsis ( ω ). El último que cierra el conjunto de elementos orbitales es la inclinación ( i ), que mide el ángulo formado entre los dos planos, contenido entre 0 y 180 grados. En la Figura 2.1 puede comprobarse el sentido positivo de cada uno de estos tres ángulos. Para resumir, queda definida en el espacio cualquier órbita no degenerada mediante los siguientes parámetros4: •p: Parámetro de la cónica. •e: Excentricidad. •θ: Anomalía Verdadera. •Ω: RAAN. •ω: Argumento de Periapsis. •i: Inclinación. Diferenciación entre distintas cónicas Se ha demostrado al comienzo del Apartado 2.1.1 la primera Ley de Kepler, que enuncia que la trayectoria seguida por un cuerpo sumergido en la acción gravitatoria de otro es una cónica; así como se ha expuesto en función de e la cónica en cuestión. En este apartado, se realiza un breve repaso por cada una de estas cónicas, haciendo énfasis en lo que presenta mayor utilidad para el objetivo del trabajo: la relación del parámetro de la cónica ( p ) con la curva en cuestión según su tipo. En el primero y más sencillo de los casos, la circunferencia ( e=0 ), la relación es clara y directa: r=p 1+ecosθe=0 =p,(2.15) con lo que para este caso, el parámetro pno es más que el radio de la circunferencia. Figura 2.2 Esquema de órbita elíptica en plano perifocal, extraída de [12]. 4 Los elementos orbitales no son únicos, dependiendo del autor puede encontrarse alguno de estos intercambiados. Para este trabajo, estos han sido los escogidos.
16 Capítulo 2. Modelo e hipótesis Las Leyes Horarias son las siguientes6: •Elipses: 1.x | tan(θ/2) = q1+e 1−etan(E/2). 2.x | M=E−esinE. 3.x | M=n∆t. •Hipérbolas: 1.y | tan(θ/2) = q1+e e−1tanh(H/2). 2.y | N=esinhH−H. 3.y | N=n∆t. •Parábolas: 1.z | 2B=3tanθ 2+tan3θ 2. 2.z | B=3qµ p3∆t. 3.z | z=3 pB+√B2+1. 4.z | tanθ 2=z−1 z. El procedimiento de resolución es el siguiente: 1- Dada una posición que no corresponda con θ=0 o θ=180º ( θ∞ en caso hiperbólico, véase la Figura 2.4), en cuyos casos los tiempos desde periapsis son conocidos (0 en el primer caso, T/2 en el segundo si es una elipse o infinito si es una parábola o hipérbola), se resuelven las ecuaciones en el orden (1.x)-(2.x)-(3.x) (respectivamente, (1.y)-(2.y)-(3.y) en el caso hiperbólico y (1.z)-(2.z) en el parabólico) para determinar el tiempo que corresponde hasta ese punto desde periapsis (el punto puede ser posterior o anterior a este cruce, los tiempos de vuelo son simétricos). 2- Se calcula el tiempo a incluir en las ecuaciones para determinar la posición final, en el que pueden darse las siguientes casuísticas: – El ángulo de la posición inicial se encuentra entre 0 º y 180 º , en cuyo caso el tiempo de cálculo será la suma del calculado en el punto 1- más el tiempo de propagación. Si se está en el caso elíptico y el tiempo calculado supera el semi-período, es recomendable aplicar propiedades de simetría y calcular el ángulo correspondiente a un tiempo desde periápsis de T−∆t , siendo el ángulo de la posición final 360 º menos el valor que den las ecuaciones. – El ángulo de la posición inicial es superior a 180 º en una elipse o negativo en los otros dos casos, y el tiempo desde dicho punto hasta periapsis es menor que el tiempo de propagación. El tiempo a usar será entonces la diferencia entre el de propagación y el calculado en 1-, y el ángulo obtenido será directamente el de la posición final (si se diera el caso en el que ∆t sea mayor al semi-período de una órbita elíptica, puede seguirse la misma estrategia que en el caso anterior). – El ángulo de la posición inicial es superior a 180 º en una elipse o negativo en los otros dos casos, y el tiempo desde dicho punto hasta periapsis es mayor que el tiempo de propagación. En este caso el ∆t a usar en las ecuaciones será el tiempo calculado en 1- menos el de propagación, y el ángulo final será 360 º menos el obtenido (o simplemente menos el ángulo en hipérbolas y parábolas). 6(2.z) y (2.x) son denominadas, respectivamente, Ecuación de Barker yEcuación de Kepler.
2.1 Conceptos fundamentales de la Mecánica Orbital 17 3- Con ∆t se resuelven las ecuaciones en la dirección contraria ((2.z)-(3.z)-(4.z) para parábola, siendo indiferente la solución que se elija en z sobre el valor final de θ ). Cabe destacar que (2.x), la Ecuación de Kepler, es una ecuación trascendental, sin solución analítica. Esta ecuación y (2.y) deben resolverse utilizando algún método numérico. El problema puede enfocarse también en el sentido contrario de resolución, dadas dos posiciones en la órbita determinar el tiempo entre ellas, resoluble con mayor facilidad al no tener que resolver las ecuaciones trascendentales, pero carece de utilidad tal y como se resuelve el problema planteado en este proyecto. 2.1.3 Problema de Lambert Se ha expuesto cómo determinar una órbita dados el vector de posición y velocidad en un instante, o cómo propagar en el tiempo la posición dentro de una órbita. El Problema de Lambert presenta el objetivo de determinar una órbita a partir de dos posiciones conocidas , y el tiempo que se tarda en ir de una a otra siguiendo la órbita que se pretende encontrar. Las ecuaciones que aparecen para resolverlo no son triviales, y trabajos de investigación completos se dedican a encontrar la manera óptima de hacerlo. Existen multitud de algoritmos distintos, cada uno con sus ventajas y desventajas. En este trabajo se ha utilizado el procedimiento presentado por Bate, Mueller and White (1971) and Bond and Allman (1996) [ 2 ]. Más detalles sobre el algoritmo pueden apreciarse en el apéndice B del documento, este apartado se dedica para plantear el problema. Figura 2.8 Esquema del problema de Lambert, extraída de [2]. En la Figura 2.8 pueden apreciarse algunas de las variables del problema. En un Sistema de Referencia similar al descrito en este documento en la Figura 2.1, se aprecian el vector de posición inicial r1 y final r2 , así como la Anomalía Verdadera (el incremento) medida en el sentido de avance del movimiento. Aparece también un nuevo parámetro c que representa la cuerda del arco de trayectoria.
18 Capítulo 2. Modelo e hipótesis A estas variables se añade el tiempo de vuelo entre ambos puntos y el número de revoluciones permitidas (parámetro que no se contempla en este trabajo), y es posible determinar la órbita que cumple las condiciones. Para elegir correctamente ∆θ, que se calcula como: cos∆θ=r1·r2 r1r2 ,(2.34) es suficiente con calcular el producto vectorial de los dos vectores ( r1∧r2 ), tal que si la componente en z del vector resultado es negativa, indicará que el ángulo correcto que se corresponde con el movimiento es 360 º menos la primera solución del coseno (se parte de r1 ). Hay que hacer también la distinción de la situación en la que la órbita fuera retrógrada, en cuyo caso el criterio es justo el contrario. Siguiendo el procedimiento expuesto en [ 8 ] se pueden obtener las ecuaciones clásicas del Problema de Lambert. Haciendo uso de las leyes horarias, relaciones geométricas entre los vectores posición y E o H según el caso, puede llegarse a los siguientes sistemas de ecuaciones que permiten resolver el problema: 1 Caso elíptico: √16µ (r1+r2+c)1/32∆tsin3(α/2) = α−β−sinα+sinβ; 0 ≤α≤2π,(2.35) sin(β/2) = qsin(α/2);−π≤β≤π,(2.36) q=±√r1r2 (r1+r2+c)/2cos(∆θ/2).(2.37) Se debe elegir la solución positiva de q (conocido con los datos de entrada al problema) si ∆θ está entre 0 y π , y la negativa en caso contrario. Una vez obtenidos los valores de α y β es posible determinar el valor del semieje mayor de la elipse como: c 2a=sin((α+β)/2)sin((α−β)/2).(2.38) Conocido a , e puede obtenerse haciendo uso de (2.1) que está expresada en función de p(a,e) , θ y ecreando un sistema de ecuaciones con r1yr2, sabiendo ∆θ. De forma análoga pueden obtenerse las ecuaciones que permitirían obtener la cónica en caso de ser hiperbólica o parabólica: 2 Caso hiperbólico/parabólico √16µ (r1+r2+c)1/32∆t=−(sinh2(γ/2))−3/2[γ−δ−(sinhγ−sinhδ)],(2.39) sinh(δ/2) = qsinh(γ/2),(2.40) q=±√r1r2 (r1+r2+c)/2cos(∆θ/2),(2.41)
2.2 Peculiaridades de las transferencias interplanetarias 19 donde γ y δ juegan papeles similares (aunque con particularidades) a los de α y β en el caso elíptico. El caso de γ=0 se corresponde con la parábola, que requiere un especial tratamiento de las ecuaciones y desarrollos en serie para evitar errores numéricos elevados en su resolución. En la búsqueda de una trayectoria óptima (objetivo) sólo aparecerán casos hiperbólicos cuando los tiempos de vuelo exigidos sean bajos, por ello no se hace una mayor profundización en el uso de γ y δ . Esto no quiere decir que el caso no se haya tenido en cuenta, y se puede ver cómo se trabaja con él y cómo se obtiene la cónica en el apéndice B. 2.2 Peculiaridades de las transferencias interplanetarias El problema de estudio se centra en la transferencia entre dos planetas del Sistema Solar. Por ello, es lógico escoger para trabajar el Sistema de Referencia Heliocéntrico Eclíptico , que tiene por características: • Origen del sistema de referencia en el Sol (muy aproximadamente el foco de todas las órbitas planetarias del Sistema Solar). •Eje x apuntando hacia el Primer Punto de Aries à. •Plano de referencia: plano de la eclíptica (en el que está contenido la órbita terrestre). •Eje z tal que las órbitas de todos los planetas son directas. Figura 2.9 Sistema de Referencia Heliocéntrico Eclíptico. Las órbitas de los planetas del Sistema Solar han sido estudiadas al detalle y con mediciones desde los tiempos de Kepler, por ello en lugar de realizar una hipótesis demasiado estricta como sería suponerlas circulares, para aportar realismo al proyecto se ha hecho uso de las Efemérides Planetarias. 2.2.1 Efemérides Planetarias Las efemérides escogidas [ 11 ] están basadas en el Sistema de Referencia J-2000, y cuentan con gran precisión para los años 1800-2050 d.C.; y con menor precisión (lo cual no quiere decir que
20 Capítulo 2. Modelo e hipótesis sean erróneas las medidas, están de hecho corregidas con términos extras) para los años desde 3000 a.C. hasta 3000 d.C. Estas efemérides trabajan con un modelo de perturbaciones en el que los elementos Keplerianos (los expuestos en esta memoria) se actualizan según el tiempo correspondiente siguiendo un modelo lineal de la forma: c=c0+˙ c·S,(2.42) donde c es un vector que recoge los seis elementos de Kepler de la órbita de un planeta, c0 los valores de estos elementos en el modelo J-2000, y S los siglos pasados desde el año 2000, calculados a partir de Días Julianos7tal que: S=t−2451545.0 36525 .(2.43) Una aclaración importante es que dada la hipótesis de problema plano, se ha optado por imponer la inclinación de las órbitas de todos los planetas (y sus variaciones) nulas. Existen otros métodos de llevar a cabo esta hipótesis, como proyectar las órbitas en el plano de la eclíptica. Dados los bajos valores en inclinación de las órbitas de los planetas que se escogen para las simulaciones ambas aproximaciones son muy similares, este método sería sin embargo poco realista si se fuera a analizar una misión a Plutón (aunque no sea un planeta), cuya órbita se sitúa en torno a los 17 º de inclinación, si bien igualmente válido para un análisis de misión preliminar. 2.2.2 Esfera de Influencia Un concepto importante de cara a trabajar en el Sistema Solar bajo la hipótesis de los dos cuerpos, es el de Esfera de Influencia. En los programas desarrollados no se trabaja explícitamente con ella, pero sí que sirve para justificar una de las hipótesis en la que se basan, el ajuste de cónicas. La Esfera de Influencia se puede definir como el espacio geométrico asociado a un cuerpo masivo en el que en un problema gravitatorio con dos cuerpos masivos y uno no masivo, la influencia de uno de los dos (masivos) puede considerarse una perturbación frente a la del otro (al que está asociado la esfera). El radio de la Esfera de Influencia respecto al Sol de un planeta cualquiera P, puede calcularse como: RE.I.−P=DÀ−Pµp µÀ2/5 ,(2.44) siendo DÀ−P la distancia del Sol al planeta, y µi=G·mi . Si se evalúa este radio para los distintos planetas del Sistema Solar, se obtiene la siguiente tabla (extraída de [12]): 7 Los Días Julianos son una medida de tiempo definida por Joseph Salinger en 1582 como una cuenta de días a partir de las 12:00 UT del 1 de Enero del año 4713 aC; que busca evitar ambigüedades como años bisiestos en escalas de tiempo astronómicas, permitiendo unificar calendarios [12]. Un día Juliano consta de 24 horas.
2.3 Variables de optimización de la misión 21 2.3 Variables de optimización de la misión En un problema de optimización es crucial tanto escoger las variables que entran en juego, como el tratamiento que se les da, esto último es especialmente importante en cuanto a al convergencia de resultados se refiere. Según los autores de la bibliografía consultada, tanto variables como tratamientos presentan mucha variabilidad. Para este trabajo se ha elegido introducir las maniobras de espacio profundo (modeladas como impulsos instantáneos), los tiempos asociados a estas, y la velocidad y fecha de salida. En el caso en el que el número de maniobras intermedias fuera nulo, las única variables serían el tiempo total de la misión y la fecha de inicio, pues la velocidad de salida quedaría fijada por las posiciones del planeta de salida, el de llegada y el tiempo de vuelo, que conforman un problema de Lambert. Este es un problema típico de transferencia interplanetaria al que están asociadas las denominadas Porkchop Plots. Estas gráficas y el problema son evaluados para verificar el funcionamiento de los programas herramientas y el algoritmo genético al comienzo del Capítulo 5 de Simulaciones y Resultados. El tratamiento de las variables se reserva al Capítulo 3 en la que se desarrolla el algoritmo genético, donde se puede entender mejor el por qué de las decisiones tomadas. 2.3.1 Maniobras de Espacio Profundo Una Maniobra de Espacio Profundo (Deep Space Maneuver) consiste en, en un punto a elegir de la trayectoria interplanetaria de un satélite, realizar una maniobra de carácter instantáneo que cambia la velocidad del cuerpo, en módulo, dirección, o ambos, sin que este cambie su posición. De esta forma, con el impulso el satélite modifica su órbita sin modificar su posición en el instante en que se produce la maniobra. Figura 2.10 Esquema de Maniobra de Espacio Profundo. Su utilización es tal y como se muestra en la Figura 2.10: tras un tiempo determinado en la órbita anterior, se aplica en la posición final (para lo que son necesarias las leyes horarias) el impulso; generándose así una nueva trayectoria que se recorre durante el tiempo hasta la siguiente maniobra, o hasta llegar al planeta destino.
22 Capítulo 2. Modelo e hipótesis Para evitar que tras el último tiempo después de la última maniobra no se alcanzara el planeta destino, la última maniobra de espacio profundo (penúltima maniobra del problema, siendo la última la de aparcamiento en el planeta destino) no es una variable, sino que se calcula a partir del Problema de Lambert tal que en la posición en la que se aplica se consiga la velocidad que cierra la misión uniendo este punto y el planeta destino (en su posición en la fecha de llegada) en el tiempo restante. Este tipo de maniobras se utilizan en la industria y han permitido el éxito de numerosas misiones, no sólo a nivel de optimización, sino como método de corrección de las órbitas. En artículos como [ 10 ] o [ 1 ] puede verse que es una variable típica de optimización, que se puede intercalar además con Maniobras Asistidas por Gravedad mejorando su factibilidad, o incluso utilizarse para mejorar el Rendezvous en una órbita. 2.3.2 Velocidad y fecha de salida La fecha de salida determina a partir de las efemérides planetarias la velocidad y posición del planeta de salida, y junto a la suma de los tiempos (la fecha final) la posición y velocidad del planeta de destino. Figura 2.11 Velocidad de salida del planeta de salida. La velocidad de salida, directamente relacionada con uno de los impulsos de la misión, es otra de las variables del problema. Esta velocidad se corresponde, tal y como puede observarse en la Figura 2.11 con la velocidad en el infinito (de salida de la Esfera de Influencia) de una hipérbola que se genera a partir de una órbita de aparcamiento con un impulso tangente. La hipérbola que en en el esquema se muestra con salida paralela a la velocidad del planeta, no tiene por qué tener esa forma, y es parte de la variable el ángulo con el que escapa relativo a la velocidad del planeta.
2.4 Hipótesis realizadas 23 2.4 Hipótesis realizadas Las hipótesis que se van a tomar han sido todas mencionadas en los desarrollos anteriores. En este apartado se van a explicar en mayor detalle, si bien ya se han visto todas sus consecuencias. 2.4.1 Problema de los dos cuerpos Una de las hipótesis fundamentales aplicadas es la de suponer un problema de dos cuerpos, es decir, considerar que en el problema de estudio sobre el satélite sólo actúa un campo gravitatorio al mismo tiempo, suponiendo el efecto del resto de cuerpos del Sistema Solar así como la presión de radiación solar una perturbación en este caso despreciable (aunque sí se han tenido en cuenta en las órbitas planetarias). Esta hipótesis tiene cabida dentro del estudio de carácter preliminar, pues lo que se busca es un punto óptimo de partida. Como ya se ha explicado con anterioridad, esta y el resto de hipótesis deberán ir relajándose una vez encontrada la solución, realizando modificaciones sobre la misma para adaptarla al problema real. Se recuerda que lo que en este trabajo se aborda es el primero de los pasos de optimización. Refuerza esta hipótesis el concepto de Esfera de Influencia y la tabla con resultados de la misma expuesta en su correspondiente apartado. En ella puede observarse el radio de dicha esfera es mucho mayor que el del planeta evaluado, permitiendo considerar un espacio teóricamente infinito a su alrededor (es decir, se consideran los límites de la Esfera de Influencia como el infinito); y a su vez el tamaño de este espacio es mucho menor que la distancia que existe entre el planeta en cuestión y el Sol, provocando que durante el trayecto interplanetario sea este astro el que domina el problema. La hipótesis del problema de los dos cuerpos no queda lejos de la realidad, permitió a Kepler enunciar sus leyes y es ampliamente usada tanto en diseño como control de múltiples misiones, realizando a posteriori correcciones como las utilizadas en las efemérides planetarias aplicadas en este trabajo. Cuerpo Masivo Dentro de la hipótesis anterior, se realiza además la aproximación (perfectamente válida) en la que el planeta o astro es mucho más masivo que el satélite de estudio. De esta forma, el parámetro µ=G(mp+ms)≈G·mp, y el astro en cuestión se sitúa en el foco de la cónica. 2.4.2 Ajuste de cónicas El ajuste de cónicas se apoya en el problema de dos cuerpos, y consiste en dividir la trayectoria interplanetaria en partes, actuando en cada una de las partes sólo dos cuerpos. Las distintas partes están dominadas por planetas (de salida o de llegada), o por el Sol. El método consiste en cambiar de Sistema de Referencia al cambiar entre partes del problema, para lo cual, se supone que los Sistemas de Referencia asociados a los planetas son Inerciales, que si bien no es cierto, los errores cometidos por aceleraciones son muy bajos, y quedan dentro del error que presenta el problema de por sí por existir el resto de hipótesis. Para realizar este cambio se hace uso del concepto de Esfera de Influencia, y tal y como se muestra en la Figura 2.11, la velocidad en el infinito de la órbita de escape del planeta se sumará a la del propio planeta para realizar el cambio de Sistema de Referencia. Esto implica, que tanto llegadas como salidas serían siempre mediante hipérbolas o parábolas, como es lógico si se tiene en cuenta
24 Capítulo 2. Modelo e hipótesis que para escapar de la Esfera de Influencia debe llegarse al infinito. La parábola sin embargo, es una situación peculiar dado que implicaría que el satélite "queda atrapado" con velocidad nula en el límite de la esfera. Durante el trayecto interplanetario, donde se hace uso del Sistema de Referencia Heliocéntrico, las Esferas de Influencia se consideran de carácter puntual ocupando la posición del planeta, de tal forma que cuando el satélite alcanza una de estas posiciones, se vuelve a realizar el cambio de Sistema de Referencia. En este caso el cambio en la velocidad será tal que la velocidad de llegada (por el infinito, el límite de la esfera) al planeta será igual a la del satélite en el Sistema de Referencia Heliocéntrico menos la del planeta en el mismo sistema. En cuanto a la posición del vector velocidad dentro de la Esfera de Influencia, si bien queda determinada en qué mitad de la esfera se encuentra por el sentido de la misma, se supone elegible la posición de forma que el periapsis de la hipérbola de llegada es un parámetro a escoger. Esto en la realidad sería una de las pequeñas correcciones finales que habría que realizar para conseguir que se cumpla la órbita escogida. 2.4.3 Instantaneidad de los impulsos Una de las hipótesis tomadas está referida a una variable de optimización, y al problema en general. Consiste en considerar los impulsos realizados sobre el satélite como instantáneos, repercutiendo sobre el mismo en un cambio en su velocidad que no conlleva tiempo. No se considera en este trabajo otros tipos de propulsión como pudiera ser la propulsión eléctrica o cualquier otro tipo de propulsión continua, que conllevarían un tratamiento completamente distinto del problema. 2.4.4 Partida desde órbita de aparcamiento El lanzamiento conforma otro problema a optimizar, pero queda desvinculado del tratado en este trabajo en cuanto a variables de optimización (lanzamientos multietapas, hora de lanzamiento, etc.). Para este trabajo se asume que el lanzamiento ya ha sido realizado y el satélite se encuentra en una órbita de aparcamiento inicial, que se supondrá baja para simular una continuación del lanzamiento. 2.4.5 Problema plano La última de las hipótesis es una no tan necesaria, sino simplificativa. El hecho de reducir el espacio a bidimensional reduce considerablemente el número de variables a optimizar, y teniendo en cuenta la relativa planitud del Sistema Solar y el carácter preliminar del diseño resulta adecuado su uso. Esta hipótesis si bien es una ayuda, presenta también desventajas. Por ejemplo, se pierde potencial de optimización en cuanto al uso de Maniobras de Espacio Profundo, especialmente útiles en cambios de plano. Pero por otra parte, permite determinar de manera más sencilla la validez de cálculos realizados dada la existencia de soluciones teóricas óptimas, como la maniobra de transferencia tipo Hohmann. Por último, destacar que la facilitación se aplica a la resolución del algoritmo genético al trabajar con menos variables, y que la mayoría de los programas creados cuentan con la opción de trabajo en tres dimensiones, dejando así preparado el terreno para un posible trabajo futuro en el que se liberara esta hipótesis; una vez comprobada la eficacia del algoritmo genético.
2.5 Ampliación: Problema tridimensional 25 2.5 Ampliación: Problema tridimensional La consideración del problema plano como se ha comentado no es tanto una necesidad como una simplificación del problema para ahorrar en tiempos de computación. Los resultados obtenidos con esta hipótesis son en muchas ocasiones suficientes para esta primera etapa de diseño de misiones, y en otras muy próximos al caso real tridimensional en cuanto al Sistema Solar se refiere, debida a la baja diferencia de inclinación entre las órbitas planetarias. Por otra parte, la hipótesis bidimensional puede ocultar una importante ventaja que ofrecen las Maniobras de Espacio Profundo respecto a las misiones sin impulsos intermedios: los cambios de plano. Por este motivo y por verificar el funcionamiento de todos los programas elaborados (incluyendo el algoritmo genético) se decide ampliar el alcance del proyecto e incluir también una simulación de un problema tridimensional. En lo relativo a este capítulo de introducción teórica, es necesario realizar una serie de aclaraciones respecto a lo expuesto hasta esta sección. Se va a tratar de completar estos conceptos básicos de la Mecánica Orbital añadiendo aquellas excepciones y casuísticas no consideradas con anterioridad por la hipótesis de problema plano. Siguiendo el orden de exposición de contenidos, la primera aclaración necesaria es relativa a los casos degenerados de los elementos orbitales clásicos . Al considerar la posibilidad de i distinto de 0 o π , surge una nueva excepción correspondiente a una órbita circular de inclinación no nula, es decir e=0con i∈(0,π). En este caso θyωno están bien definidas al no existir línea de ápsides. Se define para solventar este problema el elemento orbital argumento de la latitud: u=θ+ω;(2.45) ángulo referido a la línea de nodos, con sentido positivo en el sentido del movimiento. Su cálculo es simple a partir del vector unitario de la línea de nodos: cosu=r·n r,(2.46) escogiéndose la segunda solución del coseno cuando la componente en zde res negativa. La segunda diferencia respecto al problema plano reside en las matrices de giro para pasar de elementos perifocales al Sistema de Referencia usado, como se expone en (2.33). En el caso más general la matriz de giro toma la forma: C= cosωcos Ω−sinωsinΩcosisinΩcosω+cosΩsinωcosisin ωsini −cosΩsin ω−sinΩcosωcosi−sinΩsinω+cosΩcosωcosicos ωcosi sinΩsin i−cosωsinicosi .(2.47) El resto de teoría y desarrollo expuesto es consistente con el problema tridimensional, pero sí aparece diferente tratamiento en las efemérides planetarias, al dejar de tomarse i=0 ; y en las Maniobras de Espacio Profundo, pasando estas a estar formadas por tres componentes en el problema en 3D. La velocidad de salida tiene también tres componentes en este escenario, pero se considera como hipótesis que la órbita de aparcamiento inicial está situada en el plano de la hipérbola de escape determinada en la solución alcanzada.
32 Capítulo 3. Algoritmo genético 3.2.2 Partes de un algoritmo genético Entendida la base del funcionamiento de un algoritmo genético, se explica en este apartado cómo se traduce esa teoría en el cálculo computacional. Se van a introducir cada una de las partes que componen el algoritmo, los denominados como operadores genéticos, qué grados de libertad existen en cada uno de ellos y su implicación, y cuáles de entre las distintas posibilidades de funcionamiento se adapta mejor al problema a resolver. Toda la codificación, realizada en MATLAB , puede verse en el Apéndice A. Generación de la población inicial La opción elegida para la generación de los individuos que conforman la primera población es realizar este proceso de manera puramente aleatoria, en términos computacionales. Esto se traduce en no fijar una semilla de generación para el proceso de asignación de números aleatorios. La semilla es una variable de entrada que determina en cualquier algoritmo de generación de números aleatorios la forma de crearlos. En programación, la tarea de proveer números aleatorios de la forma en la que lo haría un dado no es trivial, en realidad estos algoritmos no son completamente aleatorios y por ello, con una misma semilla cada vez que se comience a usar el algoritmo la secuencia de números mostrada es la misma. Para algunos problemas y para el proceso de creación del algoritmo puede ser interesante fijar todas las semillas de todos los procesos que conlleven aleatoriedad, privando al método de la eficacia que conlleva dejarlos libre pero siendo útil para analizar si la ejecución del mismo está siendo correcta. El algoritmo que se utiliza para obtener todos los resultados expuestos en el trabajo no tiene fijada ninguna de estas semillas. Para determinar estos individuos, se utiliza la orden rand de MATLAB sin fijar la semilla, asignando a cada gen de cada individuo de forma aleatoria un valor comprendido entre su máximo y mínimo. Crear un algoritmo genético que ofrezca siempre el mismo resultado, aunque este se haya conseguido en una primera iteración siguiendo el camino de la metaheurística, no tiene sentido. La gran exploración que ofrecen estos métodos se ve limitada si no deja libertad en la inicialización del problema. Un algoritmo genético no debe lanzarse una única vez, debe hacerse en varias ocasiones en pos de determinar la solución más cercana al óptimo global. Nunca es posible asegurar que se ha encontrado, pero si es posible comparar entre varias iteraciones cuál está más cerca. Selección El Operador de Selección tiene como misión elegir, tras realizar el resto de operaciones genéticas y tenerse un número de individuos Y>X , siendo X el número de individuos que conforman una población (constante), los (Y−X)individuos que serán eliminados. La selección de individuos que sobreviven y pasan a formar parte de la siguiente población es quizás la parte del algoritmo en la que la decisión de cómo afrontarla produzca mayor variedad en el resultado. Constituye la capacidad de supervivencia en sí misma, y participa en gran medida en la presión sobre la población[ 1 ]. Este es un concepto que está directamente relacionado con la exploración y explotación del método. Con unas condiciones de supervivencia estrictas, por ejemplo, si sobreviven siempre los X mejores individuos (tras los procesos de cruce y mutación, este procedimiento es conocido como Selección por Truncado) la presión será alta y con ello la explotación, pero se perderá capacidad de exploración llevando el método a una convergencia rápida a un óptimo local. Lo contrario, una presión demasiado baja como pudiera ser la selección aleatoria de individuos que sobreviven, carecería de capacidad de explotación.
3.2 Base y funcionamiento de los algoritmos genéticos, Selección Natural 33 El objetivo es conseguir un método capaz de aplicar una presión baja al principio favoreciendo la exploración, y alta en sus últimas iteraciones para explotar la mejor zona encontrada. Para lograr este objetivo existen distintos procedimientos [9], algunos de ellos son: •Fitness Proportionate Selection (FPS) : A veces referida como Selección de Ruleta, la probabilidad de supervivencia de un individuo es proporcional a su valor en la función objetivo. Esta probabilidad viene determinada en el caso de maximización por la expresión: pi=FVi ∑XFVi ,(3.1) que cambiaría a (1−pi) en minimización, siendo FVi el valor de la función objetivo para el individuo i . Se trata de un procedimiento con presión alta, pero que permite a cualquier solución la oportunidad de sobrevivir, pues nunca se alcanza una probabilidad estricta de 0. •Stochastic Universal Sampling (SUS) : Consiste en una variación de FPS que logra disminuir la presión. El procedimiento consiste en elegir un número n dentro de Y , ordenar las soluciones según su valor en la función objetivo, y proceder a determinar cuales sobreviven utilizando probabilidades similares a las del método anterior, pero de forma que se asegura la elección de al menos un individuo dentro de cada rango (orden de magnitud) de los valores de función objetivo. De esta forma se asegura que existe siempre posibilidad de exploración. F.Cacciatore yC.Toglia presentan en [ 1 ] alternativas como Ranking Selection ySigma Scaling; pero están de acuerdo con Sean Luke, autor de [ 9 ], en que el procedimiento más adecuado para un algoritmo genético (y en concreto adecuado para el tratamiento de problemas de Mecánica Orbital según se expone en [1]) es el conocido como Selección por Torneo. La Selección por Torneo consiste como su propio nombre indica en realizar competiciones entre un número t de individuos aleatorios del grupo Y tal que el mejor de ellos pasa a formar parte de la siguiente población y deja de participar en los torneos siguientes. Se repite este torneo X veces de forma que se mantiene el número total de individuos en la población seleccionada. El parámetro t es uno de los mencionados parámetros de control sobre el algoritmo genético, encargado de controlar la presión sobre la población durante el proceso de selección. Un valor de t=1 implicaría un proceso de selección aleatoria, carente de explotación; mientras que uno superior a X consistiría en un proceso de truncado, reduciendo la exploración. Típicamente el valor de t=2logra buenos resultados, pudiendo aumentar para crear un algoritmo más selectivo. Explorando esta idea, el algoritmo que se desarrolla en este trabajo comenzará con un valor de t=2 que aumentará en iteraciones avanzadas para incrementar la presión en estas etapas y con ello la explotación. Figura 3.4 Diagrama de proceso de Selección por Torneo.
34 Capítulo 3. Algoritmo genético Cruce El Operador de Cruce es el encargado de seleccionar a dos individuos y crear uno nuevo a partir de la recombinación de sus genes (de los valores de sus variables). Para crear esta parte del algoritmo existen distintas metodologías, la mayoría de ellas enfocadas a individuos con variables binarias. Principalmente pueden distinguirse dos caminos para abordar el proceso [9]: •Recombinación por sustitución : Consiste en elegir ciertas variables de uno de los progenitores para el nuevo individuo y completar los genes de la descendencia rellenándolos con los del otro progenitor. •Recombinación lineal : Consiste en elegir cuánto se parece la descendencia a cada uno de los progenitores, siendo los genes del nuevo individuo una combinación lineal de los predecesores. Este camino no es aplicable a un problema de variables binarias. En un problema no binario (como el presente en este trabajo) cualquiera de las dos alternativas es viable así como una combinación de ellas. La naturaleza del problema abordado sin embargo provoca que el primer método conlleve una inevitable pérdida de explotación en las etapas finales del proceso. El cambio en una sola de las velocidades de maniobra, simplemente en su ángulo, puede provocar un cambio radical en las órbitas y afectar a la convergencia. En la búsqueda de obtener un método con presión selectiva que favorezca la exploración en los inicios y la explotación en el final, se ha decidido implementar el método de recombinación lineal conocido como SBR (Simulated Binary Recombination), propuesto por K. Deb en 1998 tras la continuación de sus estudios comenzados en [ 4 ]. Esta recombinación consiste en usar una variable aleatoria u comprendida entre 0 y 1 tal que: ¯ β= (2u)1 ηc+1,u≤0.5; (3.2) ¯ β=1 2(1−u)1 ηc+1 ,u>0.5; (3.3) implicando ¯ β la forma de la descendencia, tal que si xp1 y xp2 son los progenitores, y xd1 y xd2 la descendencia (en este método siempre habrá dos nuevos individuos por cada pareja de progenitores), se tiene: xd1=1 2[(1+¯ β)xp1+(1−¯ β)xp2],(3.4) xd2=1 2[(1−¯ β)xp1+(1+¯ β)xp2].(3.5) El parámetro ηc que aparece en (3.2) y (3.3) es otro de los denominados en este documento como parámetros de control, en este caso consiste en un número real no negativo que controla la probabilidad de que la descendencia sea similar a los progenitores. El valor concreto se decide en el Capítulo 4 del documento tras simulaciones variando los distintos parámetros. Lo que se consigue con SBR es el objetivo buscado con el cruce. Durante las primeras etapas donde los individuos son muy distintos entre ellos, este método de recombinación favorece la exploración del espacio completo de soluciones; mientras que a medida que el método converge y los individuos comienzan a parecerse, independientemente del valor de ηc la descendencia se parece más a los progenitores, explotando el espacio final de búsqueda. La selección de los progenitores es otro tema abierto en la creación del operador. Desde el uso de torneos para su elección hasta la elección aleatoria de los mismos, los procedimientos existentes son similares a los de la selección de población superviviente. El método que se ha elegido en este trabajo (de invención propia) es el siguiente:
3.2 Base y funcionamiento de los algoritmos genéticos, Selección Natural 35 1- Se recorre la población de individuos, ordenada de mejor a peor individuo, decidiendo con una probabilidad pc si un individuo pasa a ser progenitor o no. Cuando se consiguen dos progenitores, se produce el cruce. Esta forma de proceder favorece la explotación y convergencia promoviendo el cruce entre mejores individuos. 2- De forma independiente al paso anterior, se fuerza el cruce de forma arbitraria entre individuos cualesquiera. Esto se realiza un número bajo de veces (parámetro de control nc ) y favorece que en ningún momento cese la exploración. Figura 3.5 Diagrama de proceso de Cruce. Mutación El Operador de Mutación es una herramienta crucial para la exploración. La mutación permite, en las fases del método en la que la convergencia es alta y la exploración por el resto de operadores baja, la capacidad del algoritmo de seguir explorando con la posibilidad de encontrar mejora. Esto no quiere decir que la mutación no ayude a la explotación, hacia cuál de los dos conceptos se enfoque este operador depende de la magnitud de los cambios producidos en los genes, siendo una herramienta útil para explotar el espacio final de búsqueda si los cambios producidos son pequeños. Para controlar este operador cuando los individuos están formados por genes no binarios un buen algoritmo es el de Convolución Gaussiana propuesto en [ 9 ], con una ligera modificación. El método consiste en mutar un individuo elegido modificando genes aleatorios de forma aleatoria, pero tal que las modificaciones (siempre comprendidas entre los máximos y mínimos del gen) responden a una distribución Normal de media cero y varianza controlada σ . Ajustando el valor de σ y el de la probabilidad de que un gen mute, que se denominará pm , (nuevos parámetros de control) se consigue que el procedimiento se centre más en la explotación (menor valor de σ y pm ) o en la exploración (caso contrario). La probabilidad pm puede incluso tomar el valor unidad, en cuyo caso el individuo a mutar sufre cambios en cada uno de sus genes. Este caso es especialmente útil para conseguir que el algoritmo genético nunca pare de explorar de forma arbitraria el espacio, y por ello se va a hacer uso de él en este trabajo. La modificación realizada sobre este método es en lo relativo al uso y especificación de σ . En lugar de establecer un valor de este parámetro que afecte igual a la mutación de cualquiera de los genes, se decide que el valor de la varianza en la distribución Normal dependa del gen a modificar. Mientras la mayoría de algoritmos trabajan bien cuando todos los genes tienen el mismo orden de magnitud, este cambio permite trabajar sin necesidad de preocuparse por este aspecto (lo cual es especialmente conveniente en un caso en el que se trabaja con genes que comprenden una fracción de 0 a 1, y genes que especifican el Día Juliano, de ordenes de millares). La forma de establecer el valor de σpara cada uno de estos genes será: σ2 i=UBi−LBi σc ,(3.6)
36 Capítulo 3. Algoritmo genético donde UBi y LBi son los límites superiores e inferiores de la variable del gen respectivamente, y σc el parámetro de control cuyo efecto es inverso al de σ (no se debe confundir σc con la desviación típica, se le da este nombre por su relación con la misma: a mayor valor de σc menor valor de σ y mayor explotación, a menor valor mayor exploración por el aumento de σ). El uso del Operador de Mutación tendrá dos propósitos fundamentales: 1- De forma arbitraria cada individuo de la población Y0 (es decir con descendencia incluida) se somete a la decisión de mutar o no. Si dicha mutación sucede, se aplica el método de Convolución Gaussiana para un valor de pmyσca determinar. 2- De forma arbitraria se escoge un individuo de la población Y” (población mutada) y se le aplica una mutación con valor de pm=1. Figura 3.6 Diagrama de proceso de Mutación. Elitismo El Operador de Elitismo es una herramienta no siempre usada en algoritmos genéticos, pero de especial utilidad sobre todo en casos como el que atañe en el que pequeñas variaciones en los genes pueden provocar grandes cambios en el valor de la función objetivo. Su función es simple: preservar las E mejores soluciones de la población y evitar que estas pudieran ser eliminadas durante la selección ó sometidas a mutación. El nuevo parámetro de control E debe ser al menos siempre un orden de magnitud inferior a X (población total que sobrevive) y dependerá de él, siendo en el caso de poblaciones no muy numerosas E=2 (se recuerda que se ha decidido que al menos se puede crear una vez descendencia de la élite) hasta un máximo de E=10 en las más amplias. Dadas las características de los individuos pertenecientes a este grupo, serían necesarias ligeras modificaciones sobre las Figuras 3.4 y 3.6 pues en ellas los grupos mostrados que se someten al operador son el cómputo de la población completa tras el proceso de cruce, siendo por tanto necesario extraer de estos grupos los seleccionados como élite, y se realizarían consecuentemente X−E torneos. Se ha decidido dejar así estas figuras para contribuir a su labor explicativa en las partes correspondientes y aportarles un valor más genérico, pues como se menciona no todos los algoritmos genéticos poseen Operador de Elitismo. Este problema no aparece en la Figura 3.5, donde se permite la mezcla entre individuos normales y de élite a parte del cruce aleatorio por otro lado. En el Apartado 3.2.4 puede observarse cómo queda la integración entre las distintas partes, con la especificación de qué parte de la población participa en cada operación.
3.2 Base y funcionamiento de los algoritmos genéticos, Selección Natural 37 3.2.3 Restricciones Las restricciones juegan un papel fundamental en cualquier problema de optimización, siendo en numerosos casos una barrera difícil de sobrepasar a la hora de encontrar soluciones factibles. Elegir la manera de afrontar el problema cuando es posible puede reducir e incluso eliminar en ocasiones las restricciones del problema, escoger adecuadamente las variables de entrada y la forma de proceder puede suponer como ya se ha comentado la conversión de un problema NP-Duro en NP-Completo. Principalmente pueden catalogarse en dos grupos: •Restricciones lineales: Son las más sencillas y fáciles de abordar, establece una relación lineal entre todos o algunos genes que debe cumplirse para que la solución sea aceptable. Puede expresarse matricialmente como: Ax ≤b,(3.7) siendo A la matriz de condición que relaciona los genes del individuo x , y b el vector que contiene la condición para la relación. •Restricciones no lineales: Son complejas y conllevan operaciones no lineales, lo que impide expresarlas en forma matricial. Requieren de una función auxiliar que evalúe la condición, en ocasiones esta función puede ser el propio problema y requerirse su resolución completa para determinar la factibilidad de la solución. Este tipo de condiciones es habitual en problemas NP-Duro. Pueden ser tanto de igualdad como desigualdad. El problema de estudio en este documento se ha abordado evitando restricciones no lineales gracias al cierre del mismo haciendo uso del Problema de Lambert, quedando únicamente una restricción lineal sencilla y fácil de cumplirse relacionada con los tiempos. Su forma y más detalles pueden apreciarse en el Apartado 3.4. En cuanto a su aplicación, se comprueban las restricciones cada vez que se genera un individuo nuevo, sea cual sea el operador que lo produce. Si no se cumple la restricción, se repite la operación genética hasta que se cumpla o se sobrepase un número de intentos, en cuyo caso se cesa la operación sobre el individuo y se pasa al siguiente. 3.2.4 Relación entre las partes, procedimiento En este último apartado en lo relativo al funcionamiento y bases de los algoritmos genéticos se muestra en forma de diagrama de flujo la integración entre las distintas partes que conforman el método. Como resumen previo, los Parámetros de Control2son finalmente: •X : Número de individuos que conforman la población base. •E : Número de individuos catalogados de élite. •N : Número máximo de generaciones permitido. •t: Número de integrantes de torneos. •nc : Número de cruces aleatorios por generación. •ηc : Parámetro de control de similitud descendencia-progenitores. 2 Se recuerda que la decisión sobre sus valores se reserva al Capítulo 4 tras distintas pruebas, y que puede apreciarse mayor detalle sobre la implementación numérica del algoritmo en el Apéndice A
38 Capítulo 3. Algoritmo genético •σc : Divisor en la varianza de la distribución Normal en mutación. •pmut : Probabilidad de un individuo de mutar. •pm : Probabilidad de mutar de gen en un individuo a mutar. •pc : Probabilidad de un individuo de convertirse en progenitor. •TOL : Tolerancia de discrepancia entre el valor medio en la población en la función objetivo y el del mejor individuo. Figura 3.7 Diagrama de flujo del proceso del algoritmo genético.
3.3 Pareto Dominancia, Frente de Pareto 39 3.3 Pareto Dominancia, Frente de Pareto En ocasiones en un Problema de Optimización no se persigue un único objetivo. En la práctica, el ingeniero debe enfrentarse a un problema siendo consciente de las limitaciones existentes como pueden ser el capital disponible o el tiempo admisible de duración. Estas limitaciones pueden aparecer como Restricciones al problema, o incluso imposiciones al mismo, pero es habitual que la relación entre limitantes y objetivo no sea clara y en otros tantos casos no ser una variable del problema sobre la que se tenga control. Es en este escenario donde surge la Optimización Multiobjetivo. El propósito de la Optimización Multiobjetivo como su propio nombre indica es el de optimizar al mismo tiempo dos o más problemas que están relacionados o el mismo problema enfocado a dos propósitos distintos. Una metodología típica para abordar este tipo de optimización es el de establecer una Función Objetivo única cuyo valor es una suma ponderada de los distintos objetivos. Este procedimiento en ocasiones es suficiente para obtener una buena solución, pero presenta un problema: dos soluciones muy distintas pueden tener exactamente el mismo valor en la función objetivo. Por ejemplo, si se tuviera la Función Objetivo: J=0.3·O1+0.7·O2,(3.8) donde O1 y O2 son los objetivos de optimización, y 0.3 y 0.7 sus correspondientes pesos, entonces: O1=J−0.7·O2 0.3(3.9) constituye un conjunto de soluciones para J=cte . Esto no implica necesariamente que el conjunto de soluciones sea infinito ( O1 y O2 están relacionados), pero sí que existe la posibilidad de haber varias. No se pretende decir con esto que suponga el método de pesos ponderados un mal procedimiento para resolver una Optimización Multiobjetivo, siendo en ocasiones un camino directo y perfectamente válido. Sin embargo, en un problema sobre el que se tiene poco conocimiento entre la relación de los objetivos propuestos o de las soluciones del problema en sí mismas, este método puede ocultar una posible solución mejor, al igual que ocurre en el caso de un único objetivo con los óptimos locales y globales. Para solventar la duda de si el óptimo obtenido es el mejor o no según los objetivos, surge la denominada Pareto Dominancia y con ella los Frentes de Pareto . Un Frente de Pareto no es más que una representación en el espacio de las mejores soluciones de un problema de Optimización Multiobjetivo que se construye a partir de la Pareto Dominancia, consistente en mantener en esta representación solo las soluciones no dominadas, es decir, aquellas que son al menos en uno de los objetivos mejor al resto. Este concepto puede comprenderse mejor a partir del ejemplo expuesto en las Figuras 3.8 y 3.9. En el ejemplo se muestra un problema de Optimización MultiObjetivo de una fábrica [ 9 ], que debe enfrentarse a la decisión de cuánta maquinaria elegir para sus operaciones. Producto de esta decisión surgen la energía disponible para procesar en la fábrica y el coste de mantenerlas. Se plantea entonces la maximización de la energía disponible sujeta a la minimización de costes. Para ver el comportamiento de relación entre ambos objetivos, se decide crear un Frente de Pareto (Figura 3.9) que ayude al director de la fábrica a tomar una decisión.
40 Capítulo 3. Algoritmo genético Figura 3.8 Ejemplo de Pareto Dominancia, extraído de [9]. En la Figura 3.8 puede apreciarse la denominada Pareto Dominancia. Todas las soluciones (puntos) que están dentro del rectángulo gris son dominadas por A, es decir, A es mejor en al menos uno de los objetivos a todas ellas. Estas soluciones del recuadro por tanto no formarán parte del Frente de Pareto. Esto como se puede observar, no implica que A pase a formar parte del Frente, pues queda por someterse a juicio del resto de soluciones. Figura 3.9 Ejemplo de Frente de Pareto, extraído de [9]. Una vez se han evaluado todas las soluciones y establecido la Pareto Dominancia, se forma la denominada curva de Frente no dominado de Pareto uniendo todas aquellas soluciones que no están dominadas por ninguna otra. La forma que adopta esta curva depende del problema abordado (minimización, maximización, mixto, el problema en sí mismo...). En la Figura 3.10 se exponen algunas curvas características, pero no las únicas. Es importante tener en cuenta que los ejes del Frente de Pareto no representan variables de optimización sino objetivos, en general, no coincidentes. Las soluciones (puntos) se representan por sus valores en estos objetivos, pero dos soluciones que aparecen próximas en el Frente pueden ser muy distintas entre ellas.
3.3 Pareto Dominancia, Frente de Pareto 41 Figura 3.10 Ejemplo de distintas formas de Frente de Pareto, extraído de [9]. Si bien la creación de Frentes de Pareto resulta de gran utilidad al permitir obtener una visión global del problema y aportar criterio a la decisión, no implica que deba ser siempre el camino escogido. Explorar tal cantidad de soluciones para generar el Frente conlleva tiempos elevados. Es un procedimiento que requiere justificación, un interés concreto. Los métodos usados para generar los Frentes de Pareto son efectivamente metaheurísticos, consistentes en probar una gran cantidad de soluciones pero explorando el espacio de forma eficiente para generar una curva que cubra correctamente el espectro de objetivos. Una metodología seguida ampliamente en este campo y basada en algoritmos genéticos es la conocida como NSGA-II. 3.3.1 NSGA-II Un algoritmo para generación de Frentes de Pareto aclamado por su efectividad y rapidez de cómputo es el conocido como NSGA-II . El método fundamentalmente consiste en aplicar los procedimientos de algoritmos genéticos enfocados a explorar el espacio de soluciones añadiendo a este concepto de exploración el recorrido por los distintos valores posibles de valor objetivo en pos de crear un Frente de Pareto bien definido y cubriendo el espectro de posibilidades. En este documento se entra solo en el detalle suficiente para justificar el uso de un método alternativo, la información completa puede obtenerse en [ 3 ]. Este método surge como una mejora de NSGA (Nondominated Sorting Genetic Algorithm) que busca mejorar con respecto a este en: • Mejora en los tiempos de computación a O( MN2 ). NSGA es un método complejo que conlleva tiempos de O(MN3) donde Mes el número de objetivos y Nel tamaño de población. •Añadido de Elitismo, no existente en NSGA. •Reducción de parámetros de control para crear un método más sencillo de implementar. La característica principal que diferencia a este algoritmo del resto de Optimización Multiobjetivo es la forma en la que selecciona las soluciones que permanecen en el Frente de Pareto. El número final de soluciones que conforman esta curva está limitado a un máximo de N , siendo N similar a X en el algoritmo genético, el tamaño de población de soluciones sobre el que se trabaja. Con esta limitación, puede suceder que se encuentren gran cantidad de soluciones no dominadas en un entorno no muy amplio del espectro de posibilidades, siendo entonces el Frente generado no válido para solventar el problema por el que fueron ideados: crear una visión global del problema. NSGA-II se enfoca en solventar este problema de una manera eficiente, y lo consigue con éxito como puede comprobarse en la Figura 3.11.
48 Capítulo 3. Algoritmo genético descrito (problemas tipo NP, en concreto en el problema propuesto en el proyecto) para la mayoría de aspectos de la optimización, pero no para todos y no siempre. El primer inconveniente que surge al tratar una metodología al nivel de complejidad de la presentada es la programación del algoritmo. Existen fuentes de código a disposición de cualquier usuario, el propio MATLAB contiene el suyo propio ( @ga ), pero crear algoritmo desde cero para adaptarlo al problema ofrece mejores resultados y mayor control sobre el método. En lo relativo al control surge el segundo gran inconveniente. Los algoritmos genéticos presentan un gran número de Parámetros de Control cuya calibración es necesaria para ajustar la convergencia del problema tratado. Este es un proceso que debe realizarse al inicio y puede resultar problemático ante falta de experiencia o conocimiento sobre el problema. En el ajuste previo del algoritmo se trabaja también sobre la precisión y tolerancias. Tal y como se ha expuesto en este capítulo, encontrar el equilibrio entre exploración y explotación no es una tarea sencilla, y habitualmente para conseguir tiempos de cómputo asumibles es necesario sacrificar precisión en el óptimo alcanzado. Esto hace necesaria la utilización de un algoritmo convencional auxiliar para refinar el óptimo global alcanzado, lo que infiere en un extra de tiempo empleado para resolver el problema. Otro gran inconveniente es la falta de garantías matemáticas de que la solución encontrada sea un óptimo global. Si bien la experiencia demuestra la fiabilidad de estos métodos y su gran utilidad, no es posible asegurar si quiera, previo a la resolución de un problema, que el método vaya a converger. Por último, no hay comparativa entre los tiempos de cálculo de un método convencional con uno metaheurístico. En estos últimos la necesidad de comprobar un gran número de soluciones sumada a las operaciones que se realizan con las propias soluciones entre generaciones resulta en un elevado coste computacional que se refleja en tiempos de cálculo un orden de magnitud (o hasta dos) superior a los de métodos convencionales. Esta desventaja sumada a la falta de garantías matemáticas sobre la solución alcanzada hacen de estos métodos un procedimiento inadecuado para ciertas aplicaciones que requieran rapidez y fiabilidad (lo cual no es una preocupación en el problema abordado en este documento). Estas desventajas no deben ocultar por otro lado el hecho de que los métodos metaheurísticos en ocasiones pueden resultar el único (o al menos el mejor) procedimiento para obtener el óptimo global de un problema.
4 Análisis de procedimientos del algoritmo genético Tras programar el algoritmo genético (véase el Apéndice A) se procede en este capítulo a mostrar su funcionamiento. El primer apartado se dedica a la configuración de los parámetros de control para ajustar la convergencia del método. Todos los desarrollos del capítulo se corresponden con la siguiente configuración de misión1: •Problema bidimensional2. •Misión Tierra-Saturno. •Fechas de lanzamientos comprendidas entre el 1 de Enero de 2024 y el 1 de Enero de 2025. •Órbita de aparcamiento inicial a 100 km de altitud. •Órbita final a un radio de 349392 km (seis veces el radio de Saturno). •Dos Maniobras de Espacio Profundo. •Tiempo de Vuelo Total de 600 días. Para comprobar que el algoritmo funciona correctamente, se comparan los resultados con los obtenidos por la función @ga y@fmincon de MATLAB. 4.1 Análisis de sensibilidad de parámetros de control En este apartado se realizan pruebas variando los valores de los distintos Parámetros de Control expuestos en el Apartado 3.2.4. El objetivo es primero mostrar cómo afecta al procedimiento cada uno de ellos, y finalmente escoger los más adecuados para el problema a resolver. Para visualizar su efecto se ha decidido crear una representación típica del procedimiento de un algoritmo metaheurístico consistente en representar para cada población el valor objetivo de la mejor solución (puntos rojos) y la media del resto de población (puntos negros). De esta media se ha decidido excluir un quinto de la población correspondiente a las peores soluciones, en general estas se corresponderán con la exploración más aleatoria del algoritmo y pueden aportar valores extremadamente altos en la función (incluso infinito) poco representativos de la convergencia del algoritmo durante las fases intermedias de exploración alta. 1 Se realizan los análisis con una única configuración de misión para aportar al Capítulo la posibilidad de comparar fácilmente entre los distintos cambios. Los parámetros se están ajustando para el problema, no para la configuración de la misión. Se puede apreciar que el ajuste ha sido correcto en los distintos resultados del Capítulo 5. 2A excepción del último apartado en el que se verifican los ajustes también para el caso tridimensional. 49
50 Capítulo 4. Análisis de procedimientos del algoritmo genético Se dará por válido el procedimiento del algoritmo si se cumplen las siguientes condiciones: • La curva representada por la media tiende a la curva representada por el mejor (el método converge). • La curva negra no debe ser (en general) perfectamente definida, es decir, debe presentar picos. Esto indica que se sigue explorando el espacio de soluciones. • La curva roja no debe ser constante (en general, pudiera ser que de forma aleatoria una de las soluciones generadas arbitrariamente al inicio coincidiera con el óptimo global) y a medida que avanza el algoritmo debe (en el caso de minimización abordado) tender a un valor menor al inicial. • El valor final alcanzado de la solución óptima debe ser del órden de los obtenidos con @ga y @fmincon. • Los tiempos de computación deben ser al menos del orden de los obtenidos con @ga (para esta comprobación no se pinta la gráfica de proceso). 4.1.1 Variación de X Partiendo de unos parámetros de control inicial de la forma: [X,N,TOL,E,pc,pmut,ηc,t,σc,pm,nc] = [X,400,5,5,0.5,0.3,3,2,30,0.3,2];(4.1) decididos de forma arbitraria, se procede a variar X y comprobar los cambios. El parámetro X es controlable en @ga , por lo que es posible realizar dos comparativas en para este análisis a nivel de resultados y tiempos3. Figura 4.1 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=50. 3 El problema de @ga es, como en otros tantos códigos de uso abierto, el desconocimiento sobre su programación interna. En algunos aspectos supone una caja negra y modificar sus opciones puede no repercutir en lo esperado. Este fue uno de los principales motivos que llevó a crear un algoritmo propio. Para recalcar la importancia de ajustar los parámetros de control, en @ga el número de población es el único parámetro que se modifica.
4.1 Análisis de sensibilidad de parámetros de control 51 Figura 4.2 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=100. Pudiera parecer en estas gráficas que no se está mejorando en exceso en cuanto al óptimo se refiere, es decir, que la explotación es baja. Esto es sin embargo un problema de escala, en la Figura 4.3 puede apreciarse la mejora con el paso de generaciones. La realidad es que la explotación posible depende de la aleatoriedad de la población inicial, pudiendo por suerte caer uno de estos individuos muy próximo al óptimo global. Figura 4.3 Zoom a mejora del óptimo en procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=50. Se presentan ahora unas tablas que recogen la información más relevante para comparar. En este primer análisis de parámetro se incluye @ga y @fmincon en pos de respaldar el funcionamiento del algoritmo creado y la decisión de abordar el problema por este procedimiento. Se entrará en mayor detalle en el Apartado 4.1.7.
52 Capítulo 4. Análisis de procedimientos del algoritmo genético X=50 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 758 32.0144 @ga 1162 43.4448 @fmincon 61 40.8778 Algoritmo Genético Propio (+refinador) 756 30.8784 @ga (+refinador) 1071 30.8111 X=100 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 1427 32.0584 @ga 1831 30.7426 Algoritmo Genético Propio (+refinador) 1404 32.0014 @ga (+refinador) 3839 30.7105 Con los resultados de la tabla se observa el éxito de haber programado un algoritmo genético desde cero enfocado al problema a resolver. Los tiempos son considerablemente mejores y si bien en el refinado no existen grandes diferencias entre el propio y el que incorpora MATLAB si las hay en las soluciones sin postrefinador. En definitiva de estas primeras tablas se puede concluir que se han aprovechado bien tanto la capacidad de exploración como de explotación del algoritmo. En cuanto a los métodos convencionales como @fmincon (usado en el refinado de soluciones) se observa que, pese a los escasos tiempos de cómputo que requiere, la solución alcanzada peca de depender en exceso de la solución inicial aportada (en este estudio aleatoria, excepto para el refinado), y puede converger como se ha comentado previamente a un óptimo local. No debe sorprender por ejemplo que el tiempo total sea menor en el caso con refinado que sin él para el algoritmo propio. Los tiempos son orientativos y sirven para reflejar el orden de magnitud en el que se procesa la solución, el tiempo concreto depende de la aleatoriedad implícita en el problema. En cuanto al número de individuos por generación ( X ) se puede concluir que aumentarlo en exceso repercute notablemente en los tiempos (ha de tenerse en cuenta que se va a repetir numerosas veces el uso del algoritmo para generar el Frente de Pareto), habiendo en este caso incluso para el algoritmo propio empeorado la solución convergiendo a un óptimo local 4 (se aprecia en la proximidad de la solución del genético con el refinado, y la diferencia con el caso de menor población). Por ello se concluye que la población de 50 es válida para al menos hasta dos Maniobras de Espacio Profundo en el algoritmo creado. 4.1.2 Variación de N Como se puede prever, en caso de que la solución no se alcance por tolerancia, el parámetro N determinará el tiempo de procesamiento y la cercanía de la solución final alcanzada por el algoritmo al óptimo real, tal y como se observa en las Figuras y las tablas que le prosiguen. Se varía entre N=400 yN=200, el resto de parámetros se mantiene igual a (4.1) con X=50. 4 Se recuerda que los algoritmos metaheurísticos no pueden asegurar la convergencia hacia un óptimo global, pero tienen más probabilidad de producir esta situación que los métodos convencionales.
4.1 Análisis de sensibilidad de parámetros de control 53 Figura 4.4 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y N=100 . Compárese con la Figura 4.1. N=200 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 367 35.3995 Algoritmo Genético Propio (+refinador) 424 30.8268 Las diferencias en el caso con refinado son indistinguibles, pero igual que ha ocurrido así pudiera haber sucedido que se hubiera quedado el algoritmo genético cerca de un óptimo local. Cuanto mayor sea N mayor número de oportunidades de encontrar el global. Se debe optar por un valor de compromiso entre esta seguridad y el tiempo de computación para la generación del Frente de Pareto, aunque lo ideal sería que el algoritmo encotrara la solución por tolerancia en un tiempo razonable, no alcanzando el número de generaciones N. 4.1.3 Variación de TOL Variar la tolerancia en la solución (la cual se ha definido excluyendo el peor quinto de población para comprobar si se ha alcanzado) afecta no a la convergencia, sino a los tiempos de procesamiento e indirectamente a la explotación del método. El hecho de evitar dar por concluido el método disminuyendo la tolerancia obliga al algoritmo a seguir explorando y explotando, mejorando la solución alcanzada. Si bien se mejora en el óptimo en el caso sin refinar, en el refinado para las tolerancias manejadas, las diferencias son indistinguibles.Se varía entre TOL =1 y TOL =5 , el resto de parámetros se mantiene igual a (4.1) con X=50. TOL=1 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 781 31.3304 Algoritmo Genético Propio (+refinador) 838 30.6245
54 Capítulo 4. Análisis de procedimientos del algoritmo genético Figura 4.5 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y TOL =1 . Compárese con la Figura 4.1. Si bien ya se aprecia el efecto de la tolerancia en estas tablas, esto es debido a la relación que tiene con el aumento del parámetro t en la etapa final del método. El verdadero efecto sucede cuando el algoritmo alcanza solución por tolerancia y no por número de generaciones. Mayor detalle de este fenómeno se aprecia en el Apartado 4.1.7. 4.1.4 Variación de tyE El parámetro t (participantes en los torneos de selección) resulta ser uno de los parámetros que más afectan a la explotación y convergencia del método. Tal y como se explicó en el Apartado 3.2.2, el número de participantes se aumenta cuando el algoritmo decide que se está acercando a las generaciones finales, asegurando la convergencia, mejorando la explotación en esta última etapa, y reduciendo los tiempos (la implementación puede apreciarse en el Apéndice A). Por otra parte el número de individuos de élite E presenta un efecto similar al de t , siendo cuanto mayor su valor la explotación más intensa. Para reflejar la importancia de estos parámetros, se muestran resultados para valores iniciales de t=2 con E=5 ; y t=5 con E=7 . El resto de parámetros se mantiene igual a (4.1) con X=50. t=5, E=7 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 788 32.4609 Algoritmo Genético Propio (+refinador) 958 30.9710 Como se aprecia en la Figura 4.6, el aumento de estos parámetros lleva a una convergencia excesivamente rápida del método y reduce la exploración en gran medida. En este caso se ha alcanzado una buena solución pero se recuerda que la solución alcanzada para unos mismos parámetros no es siempre idéntica y depende de la componente aleatoria del método.
4.1 Análisis de sensibilidad de parámetros de control 55 Figura 4.6 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), t=5 y E=7 . Compárese con la Figura 4.1. 4.1.5 Variación de pc,ηc, y nc El Operador de Cruce juega un papel fundamental en el proceso. Los tres Parámetros de Control que intervienen en su funcionamiento afectan sustancialmente a la convergencia. Tras varias iteraciones, los valores de los parámetros que mejor se adaptan al problema son: pc=0.4 , ηc=3 , y nc=2 . Reducir ligeramente la probabilidad de cruce reduce los tiempos sin afectar en gran medida a la exploración. El valor de ηcse selecciona de manera que los descendientes puedan llegar a ser con relativa frecuencia igual de parecidos a ambos progenitores ( ¯ β∼0.5 ). También después de varias ejecuciones del algoritmo variando nc se ha escogido este con valor de 2, lo que implica cuatro descendientes con mayor aleatoriedad, dando buenos resultados en la exploración. Figura 4.7 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), pc=0.4 , ηc=3, y nc=2. Compárese con la Figura 4.1.
56 Capítulo 4. Análisis de procedimientos del algoritmo genético pc=0.4, ηc=3, nc=2 Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 743 31.4569 Algoritmo Genético Propio (+refinador) 742 30.6877 4.1.6 Variación de pmut,σc, y pm La influencia del Operador Mutación en la convergencia del método es crucial. No sólo determina la exploración en las etapas finales permitiendo posibles mejoras sino además la convergencia gracias a los pequeños cambios controlados por σc en los genes. Las modificaciones sobre este operador resultan ser las que mayor efecto tienen sobre el método, permitiendo cerrar el problema por haber alcanzado tolerancias (gran mejora en explotación) pero con el peligro de llevar al algoritmo a óptimos locales. Figura 4.8 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), pmut =0.15 , pm=0.3, y σc=30. Compárese con la Figura 4.1. Tiempo de computación [s] ∆VTalcanzado [km/s] Algoritmo Genético Propio 649 31.8670 Algoritmo Genético Propio (+refinador) 641 31.9618 4.1.7 Decisión final sobre Parámetros de Control Como se ha podido apreciar en todas las Figuras de los apartados anteriores, ante las distintas configuraciones de los Parámetros de Control se cumplen los criterios establecidos al inicio de la Sección 4.1, validando el procedimiento del algoritmo. Tras observar el efecto de modificar cada uno de estos parámetros se realizan pruebas extras para determinar qué combinación de los mismos ofrece resultados suficientemente cercanos al óptimo global (que se puede establecer tras todas las pruebas que se encuentra en el entorno de ∆VT=30.7 km/s) en tiempos reducidos de cara a preparar el algoritmo a enfrentarse a la creación del Frente de Pareto. Tras las distintas pruebas, se ha determinado que una buena elección de los parámetros de control es: [X,N,TOL,E,pc,pmut,ηc,t,σc,pm,nc] = [50,400,3,7,0.4,0.15,3,4,30,0.3,2].(4.2)
4.1 Análisis de sensibilidad de parámetros de control 57 Con estos parámetros se lanza varias veces el algoritmo para comprobar que los resultados no son buenos por pura aleatoriedad. Se obtienen los siguientes resultados: Algoritmo Genético Propio Tiempo de computación [s] ∆VTalcanzado [km/s] Prueba 1 611 41.4235 Prueba 2 504 31.2703 Prueba 3 549 32.6524 Prueba 4 215 32.4654 Algoritmo Genético Propio (+refinador) Tiempo de computación [s] ∆VTalcanzado [km/s] Prueba 1 374 30.7162 Prueba 2 609 30.7621 Prueba 3 505 31.3902 Prueba 4 177 30.6952 Otros algoritmos Tiempo de computación [s] ∆VTalcanzado [km/s] @ga 1162 43.4448 @fmincon 61 40.8778 @ga (+refinador) 1071 30.8111 A la vista de los resultados se podría concluir en el éxito del algoritmo creado, ofreciendo mejores valores para el problema tratado en menor tiempo que @ga sin modificar sus opciones . El uso de un refinador ( @fmincon ) para la solución alcanzada por el algoritmo genético resulta más una necesidad que una opción si se busca precisión. Si bien es cierto que las diferencias en ocasiones entre ambas soluciones son bajas, las inferencias en el tiempo de cálculo no son si quiera apreciables haciendo que merezca la pena realizar la mejora. Es respecto a estos tiempos de cómputo donde se aprecia la mayor ventaja de haber programado y ajustado al algoritmo para el problema concreto. Teniendo en cuenta que se va a ejecutar multitud de veces ( m veces, siendo m el número de posibles Maniobras, para cada Tiempo Total de Vuelo definido en el Frente de Pareto) el algoritmo, cada segundo reducido por ejecución puede suponer una gran diferencia en el tiempo total de proceso. Una vez elegida la solución del Frente que se quiere realizar, lo lógico es volver a ejecutar para ese tiempo el algoritmo reduciendo su TOL y aumentando N. Figura 4.9 Procedimiento del algoritmo con los Parámetros de Control finales.
64 Capítulo 5. Simulaciones y resultados Es el conocimiento de esta solución el que va a permitir verificar el funcionamiento de los programas. Al obtener una PorkChop plot debería obtenerse un único mínimo y una forma repetida en el tiempo. Al comprobar la órbita de transferencia correspondiente al mínimo, debería observarse la tipo Hohmann. Figura 5.1 Ejemplo (extraído) de PorkChop plot para misión Tierra-Marte, caso tridimensional. La forma de resolver este problema consiste en realizar un barrido en las fechas de salida y llegada, para cada una de ellas establecer la posición de los planetas implicados, y resolver con estos datos del Problema de Lambert. Realizando esto para una misión Tierra-Saturno se obtiene las gráfica que se observa en la Figura 5.2. Figura 5.2 Porkchop plot de misión Tierra-Saturno, representando C3T.
5.1 Verificación de programas, PorkChop 65 Se procede ahora a buscar con el algoritmo genético el óptimo para comprobar si se corresponde con el que debería según las gráficas. Se trabaja sobre ∆VT por coherencia con el resto de resultados del documento, y sobre C3T para comparar. Se usa una Función Objetivo distinta a la expuesta en el Apartado 3.4, pues en este caso las variables de optimización son la fecha de llegada y salida, y en uno de los casos la función objetivo cambia (a minimizar C3T). Sin utilizar refinador, se obtiene: Caso Fecha Salida [JD-24 ·105]Fecha Llegada [JD-24 ·105]Valor Minimización ∆VT6.0447·1046.2876·10412.8713 km/s Minimización C3T6.0447·1046.2563·104137.9973 (km/s)2 Figura 5.3 Solución minimizando ∆VT.Figura 5.4 Solución minimizando C3T. Como se aprecia, ambas soluciones son casi idénticas (el algoritmo se ha lanzado en dos ocasiones distintas, una para cada Función Objetivo). Esto sucede debido a la estrecha relación entre C3T y ∆VT , el ótpimo se encuentra en las mismas fechas, pero el valor alcanzado en él de ∆VT depende de las órbitas de aparcamiento. Se comprueba que efectivamente el algoritmo ha encontrado (incluso sin refinador, no necesario por las pocas variables) la solución óptima global en el espacio de búsqueda. Además puede observarse que dicha solución se corresponde en efecto como se había predicho con una maniobra de transferencia tipo Hohmann 3. Con esta demostración se dan por válidos el conjunto de programas que intervienen en el Problema de Optimización. No sólo se ha puesto a prueba el encargado de resolver el Problema de Lambert en la generación de la Porkchop, también se ha hecho uso de propagadores, conversores de elementos a posición y velocidad (y viceversa), conversores de Días Julianos, y otros de menor importancia para realizar la representación de la misión. 3 El tiempo de una transferencia de Hohmann es TH=π n , con n=8µÀ (r1+r2)31/2 y r1 y r2 las distancias (en este caso medias) de ambos planetas al Sol.
66 Capítulo 5. Simulaciones y resultados 5.2 Misión Tierra-Saturno El primer escenario que se plantea es el que se ha venido usando para realizar pruebas sobre el algoritmo: una misión Tierra-Saturno (bajo la hipótesis de problema 2D). Esta es una misión realizada múltiples veces por distintas agencias principalmente por intereses científicos. Para el trabajo resulta útil por involucrar a un planeta tan lejano, aportando variedad a las simulaciones que se van a realizar. La configuración escogida para la misión es la siguiente: • Fechas de lanzamientos comprendidas entre el 1 de Enero de 2024 y 1 de Enero de 2025. Dada la baja excentricidad de las órbitas planetarias implicadas y la consideración de problema plano un año terrestre resulta suficiente para explorar casi la totalidad de posibles configuraciones. • Órbita de aparcamiento inicial terrestre a una altura de 100 km. Esta altura es considerada baja para una órbita de aparcamiento por la presencia atmosférica, pero su presencia en el problema existe únicamente para obviar la fase de lanzamiento. Establecer la órbita inicial es en realidad una simplificación para establecer un punto de partida, en teoría el satélite no estaría durante demasiado tiempo en esta órbita, que sería en realidad parte del propio lanzamiento. • Órbita final a un radio de 349392 km (seis veces el radio de Saturno). Se ha escogido esta órbita final para tratar de ser realista y alojar al satélite en una órbita de destino suficientemente alejada de un planeta altamente radiactivo. • Tiempos de Vuelo para la generación del Frente de Pareto que comprenden desde un año hasta veinticinco. Se han decidido estos tiempos tras comprobar en el primer apartado del capítulo que el óptimo para esta misión sin Maniobras de Espacio Profundo se encuentra en el entorno de los 5 años, y ampliando el rango pretendiendo encontrar algunas soluciones interesantes al problema. La misión más famosa lanzada a Saturno, Cassini, tardó 7 años en completar su travesía. Figura 5.5 Frente de Pareto para una misión Tierra Saturno. Primeramente, se consideran como posibilidades 0, 1 o 2 Maniobras de Espacio Profundo. Lanzando una primera vez el algoritmo genético, tras un tiempo de procesamiento de unas 3 horas, se obtienen los resultados de la Figura 5.5 (recordar que en la generación del Frente de Pareto se hace uso del refinador, por lo que lo correcto sería decir tras lanzar el algoritmo genético más refinador). Puede apreciarse una cierta tendencia clara para las misiones sin Maniobras de Espacio Profundo.
5.2 Misión Tierra-Saturno 67 Esta tendencia tiene lógica si se piensa que el Frente para el caso sin maniobras no es más que un corte a la superficie de la que se obtiene la PorkChop plot (Figura 5.2) por un plano que la cruza para cada tiempo de vuelo por el punto de mínimo ∆VT . De hecho en esta curva se aprecia un mínimo claro en el entorno de los 2000 días de vuelo, correspondiente con la maniobra de transferencia tipo Hohmann de la Figura 5.3. Esta tendencia clara en el caso sin maniobras no se reproduce en el resto de casos, donde el incremento del número de variables que intervienen en el problema ni siquiera podría representarse en un espacio tridimensional. Ampliando la Figura 5.5 en la Figura 5.6, puede comprobarse en primer lugar que en varias ocasiones las misiones convergen a unos mismos valores del óptimo. Esto se debe, como ya podía comprobarse en las Figuras 4.12 y 4.13, a que el optimizador tiene la posibilidad de en el caso de varias maniobras llevar estas a un valor nulo (o aproximadamente nulo, podrían considerarse errores numéricos). Alcanzar esta solución sin embargo no siempre es posible por el algoritmo, de ser así no existirían misiones por encima de la de cero Maniobras de Espacio Profundo (todas convergerían a la de cero), pero no es un problema en tanto en cuanto se pueda asegurar (con la seguridad que un algoritmo genético sumado a un refinador permite asegurar) que la óptima sucede sin impulsos intermedios, perdiendo interés el resto. Figura 5.6 Zoom a Frente de Pareto para una misión Tierra Saturno. Resulta lógico que en la mayoría de situaciones, y para tiempos bajos, que la solución óptima sea una trayectoria directa sin Maniobras de Espacio Profundo. Esta situación se produce por la consideración de problema plano. Es en un caso tridimensional donde se explotan las ventajas de realizar Maniobras intermedias, permitiendo cambios de plano a menor coste (véase la Sección 5.5). Haciéndose énfasis en esta primera parte del Capítulo en el problema plano, se van a analizar algunas de las soluciones que aparecen como mejores para la misión planteada en ciertos tiempos de vuelo, determinando el por qué. También se va a buscar si alguna misión mejora en ∆VT al mínimo obtenido en el caso sin maniobras. Antes sin embargo de proceder a estos análisis, conviene recordar uno de los inconvenientes (o ventajas según se mire) de los algoritmos genéticos: su resultado cambia cada vez que se lanza. Por ello se decide repetir el Frente de Pareto para tener una comparativa y validar los resultados. Se aprovecha que la velocidad alcanzada para generarlos es adecuada, y se incluye la posibilidad de hasta cuatro Maniobras de Espacio Profundo (Figuras 5.7, 5.8 y 5.9).
68 Capítulo 5. Simulaciones y resultados Figura 5.7 Frente de Pareto para una misión Tierra Saturno. En la Figura 5.7 (tiempo de procesamiento de 10 horas y media) puede comprobarse que la forma del Frente de Pareto es similar a la obtenida en la Figura 5.5, validándose los resultados. Se aprecia de nuevo la tendencia clara de las misiones sin impulsos intermedios (Figura 5.8), y aproximadamente a partir de los 3000 días de vuelo comienzan a ser mejores aquellas con Maniobras de Espacio Profundo (en este caso esto sucede para tiempos de menor valor que en la Figura 5.5, puede comprobarse que es debido a que las misiones que mejoran en estos primeros casos son las de 3 impulsos intermedios, antes no consideradas). Figura 5.8 Tendencia de misiones sin MEP en Frente de Pareto para una misión Tierra Saturno. Se vuelve ahora a ampliar la gráfica en la Figura 5.9 para seleccionar las misiones que pueda resultar de interés analizar en detalle. Se ha decidido también unir las soluciones no dominadas. La primera solución de interés es la misión con dos Maniobras de Espacio Profundo y TV=6935 JD, que resulta ser la que menor ∆VT presenta en el Frente de Pareto: 12.5709 km/s frente a los 12.6606 km/s de la tipo Hohmann. Las diferencias son mínimas pero es interesante aún así ver qué caracteriza a esta misión.
5.2 Misión Tierra-Saturno 69 Figura 5.9 Zoom a Frente de Pareto para una misión Tierra Saturno con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025. Figura 5.10 Solución con 2 MEP para TV=6935 JD. Salida el 6 de Noviembre de 2024. La solución que se ha alcanzado y puede apreciarse 4 en la Figura 5.10 es interesante, consiste en realizar algo parecido a la transferencia de Hohmann: entra en la órbita que caracteriza a esta maniobra y el satélite espera en ella tres semiperiodos hasta encontrarse con Saturno. El hecho de que la misión sea ligeramente más económica que el mínimo encontrado para cero impulsos intermedios es por el hecho de haber fijado los tiempos, por casualidad los 6935 JD en los que se desarrolla esta misión están más próximos a tres semiperiodos de Hohmann de lo que está la 4 Se representan las órbitas implicadas a escala, con los ejes en Unidades Astronómicas. La línea continua para las órbitas planetarias indica lo recorrido por el planeta en cuestión desde la fecha de salida hasta la de llegada. La línea discontinua negra señala la trayectoria seguida por el satélite. Notese la baja excentricidad de las órbitas planetarias, pero no nula.
70 Capítulo 5. Simulaciones y resultados solución sin MEP de un semiperiodo, aproximándose por tanto más la órbita a la tipo Hohmann. Es curioso sin embargo la diferencia entre la solución considerando este número de maniobras con la directa, pero la explicación es sencilla: el Problema de Lambert programado no contempla la posibilidad de múltiples revoluciones. En la solución alcanzada los impulsos intermedios podrían considerarse nulos. Sin pretenderlo, modelar varios impulsos ha solucionado un problema ni tan siquiera considerado como es la posibilidad de permanecer en una órbita de transferencia un tiempo mayor a su período. Otra solución interesante para mirar más a fondo es la única de cuatro maniobras que no está dominada, la que se corresponde con TV=8760 JD. En este caso la misión consta una primera fase de espera en una órbita elíptica para luego cercano al periapsis de esta órbita realizar el satélite una última maniobra que cambia la trayectoria para hacerlo llegar hasta Saturno. Es realmente satisfactorio el hecho de que el algoritmo determine el momento exacto en el que debe realizarse el último impulso para que este sea más eficiente5. Figura 5.11 Solución con 4 MEP para TV=8760 JD. Salida el 3 de Enero de 2024. Esta es una solución que se repite en la práctica totalidad del Frente de Pareto para las soluciones no dominadas a partir del momento en el que la directa sin maniobras deja de ser rentable. Se realiza una espera (denominada phasing) hasta cuadrar que el tiempo restante de misión sea lo más próximo al semiperiodo de la maniobra tipo Hohmann, en este instante, estando el satélite en el periapsis de la órbita de espera, se realiza un último impulso para alcanzar al planeta de destino. Este proceso se ilustra bien con la solución de 1 MEP para TV=4015 JD representada en la Figura 5.12, realizando la Maniobra de Espacio Profundo unos 2121 días tras el lanzamiento. En la Figura 5.13, solución no dominada, se puede apreciar otro ejemplo. 5Las maniobras en periapsis producen para un mismo ∆VTun mayor cambio en el tamaño de la órbita.
5.2 Misión Tierra-Saturno 71 Figura 5.12 Solución con 1 MEP para TV=4015 JD. Salida el 3 de Agosto de 2024. Figura 5.13 Solución con 2 MEP para TV=5840 JD. Salida el 22 de Septiembre de 2024. Tras analizar en detalle las soluciones es posible interpretar correctamente el Frente de Pareto. Para un problema plano en el que ambos planetas orbitan en órbitas casi circulares, la solución óptima siempre que se pueda alcanzar es la transferencia de Hohmann. Para tiempos menores al tiempo característico de esta transferencia las soluciones no pueden ser otras que hipérbolas o parábolas (la órbita de la transferencia de Hohmann es la elipse de menor tamaño que une dos órbitas circulares), en cuyo caso sigue siendo mejor dos únicos impulsos (sin Maniobras de Espacio Profundo). A partir de este tiempo de vuelo, el óptimo oscila en torno al valor de una transferencia si fueran las órbitas perfectas circunferencias, consiguiendo esto haciendo uso de varias maniobras para realizar esperas
72 Capítulo 5. Simulaciones y resultados creando la oportunidad perfecta para enlazar con el planeta de destino en una maniobra de Hohmann. Las Maniobras de Espacio Profundo por tanto se limitan en el problema plano (al menos en la primera misión planteada) a crear el mejor escenario posible para una transferencia elíptica de un solo tramo. Esta conclusión se puede observar también el trabajo de otros autores [ 1 ], donde se busca con las maniobras optimizar un Rendezvous entre dos órbitas geocéntricas que comparten plano orbital (en este caso Maniobras de Espacio Profundo no sería un nombre adecuado). No quiere decir la conclusión a la que se ha llegado que el uso de MEP carezca de utilidad en el problema planteado existiendo la posibilidad de realizar una transferencia directa tipo Hohmann, pues en ocasiones, por uno u otro motivo no se puede lanzar en la fecha correcta que lo consigue o para el tiempo de vuelo que la establece. Además, se pueden plantear escenarios interesantes que rentabilicen una misión tan lejana como es a Saturno realizando alguna aproximación a un punto intermedio que resulte de interés y regresando a la Tierra para transmitir de manera más eficiente los datos (o de mayor peso), para continuar luego su viaje hasta su planeta destino, todo ello al mismo coste que hubiera supuesto ir directamente (cualquiera de las soluciones presentadas podría atender a este escenario). 5.2.1 Análisis de explotación en generación de Frente de Pareto (Tierra-Saturno) Para terminar con esta primera misión, a modo de comprobación de cómo de bien o mal ha funcionado el algoritmo genético en lo relativo a la explotación, y verificar una vez más la validez del ajuste de los Parámetros de Control para el problema (que ya han dado buenos frutos en cuanto a tiempo de cómputo se refiere), se representa en la gráfica de la Figura 5.14 para los distintos casos calculados la diferencia entre el óptimo encontrado por el algoritmo genético y el refinado después por @fmincon . En esta representación puede comprobarse que las diferencias son inexistentes para el caso de menor número de variables, y van aumentando de forma general a medida que incrementan en número. No ha de sorprender las diferencias de hasta 10 km/s que se pueden alcanzar en ocasiones, esto no implica que el óptimo alcanzado no sea el global (un buen ejemplo de esto es que la solución en la que existe mayor discrepancia, la de 2 MEP para TV≈8400 JD, pertenece a las soluciones óptimas del Frente de Pareto). Algo que puede parecer evidente pero aún así es importante recalcar es que el refinador (usado en el Frente de Pareto) no ha empeorado las soluciones en ningún caso. Figura 5.14 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Saturno).
5.3 Misión Tierra-Mercurio 73 5.3 Misión Tierra-Mercurio El siguiente escenario a analizar será una misión Tierra-Mercurio (bidimensional). La exploración hacia este planeta inferior es famosa por sus altos costes, analizarla en este trabajo resulta de interés al involucrar al planeta más cercano al Sol, diferenciándose de la misión anterior y permitiendo explorar más la utilidad de las Maniobras de Espacio Profundo. La configuración escogida para la misión es la siguiente: • Fechas de lanzamientos comprendidas entre el 1 de Enero de 2024 y 1 de Enero de 2025. En este caso la excentricidad de la órbita de Mercurio es mayor, pero la cercanía con la Tierra hace que un año sea suficiente para evaluar multitud de casos de configuraciones iniciales. •Órbita de aparcamiento inicial terrestre a una altura de 100 km. •Órbita final a un radio de 9758.8 km (cuatro veces el radio de Mercurio). • Tiempos de Vuelo para la generación del Frente de Pareto que comprenden desde un 30 días hasta los 7 años y medio. Los tiempos han sido escogidos así para explorar tanto tiempos que incluyan travesías hiperbólicas como aquellos muy superiores donde en una revolución sea ineficiente realizar la transferencia. El caso que se estudia ahora es más complejo a nivel computacional. La mayor cercanía de las órbitas en general al Sol hace que los tiempos de búsqueda de la solución aumenten por el mayor efecto de las variables en la misma. En un caso como este el algoritmo genético es capaz de demostrar que, ante la gran diferencia de la forma de misiones, siendo el problema en cuestión el mismo, puede encontrar un óptimo global (o al menos mejor que el que aportaría un método convencional). Tras lanzar el programa y con un tiempo de procesamiento de unas 17 horas, considerando la posibilidad de hasta cuatro Maniobras de Espacio Profundo, se obtiene el Frente de Pareto de la Figura 5.156. Figura 5.15 Frente de Pareto para una misión Tierra-Mercurio con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025. 6 Se recuerda que los resultados representados se corresponden al uso de algoritmo genético propio con posterior refinado mediante el uso de @fmincon.
80 Capítulo 5. Simulaciones y resultados La peculiar forma de esta curva se explica si se analiza el Porkchop plot del problema plano de la misión Tierra-Marte entre las fechas escogidas. En la Figura 5.25 se expone este resultado, donde se puede apreciar que la forma que adopta es más parecida a la que podría obtenerse al analizarse el problema tridimensional, casi apreciándose dos tipos de mínimo distintos. Figura 5.25 Porkchop plot del C3T para una misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025. Figura 5.26 Solución con 0 MEP para TV=307 JD. Salida el 3 de Octubre de 2024. Como de costumbre se pasa ahora a analizar aquellas soluciones no dominadas más características del Frente de Pareto. Destaca que la solución que logra un mínimo ∆VT (de valor ∆VT=5.4039 km/s) se corresponde con una sin Maniobras de Espacio Profundo para TV=307 días, próximo al tiempo de vuelo de la misión Curiosity (diferencia principalmente causada por la consideración de problema plano y fechas de lanzamiento posibles). Como era de esperar, en la Figura 5.26 se
5.4 Misión Tierra-Marte 81 aprecia que esta solución se corresponde con una maniobra similar a la tipo Hohmann enlazando con Marte en el apohelio de su órbita (la órbita de Marte tiene una excentricidad no despreciable). A partir de este mínimo, en torno los 400 días de tiempo de vuelo, las misiones sin impulsos intermedios dejan de ser rentables. En la Figura 5.27 puede comprobarse que dominan las soluciones con una única maniobra intermedia en la práctica totalidad del Frente, creando una tendencia en las soluciones no dominadas bien definida y algo distinta a la de los casos anteriores, no alcanzándose en este ese mínimo en torno al cual se oscilaba en las misiones Tierra-Saturno o Tierra-Mercurio. Superados los 600 días de vuelo, las soluciones con mayor número de maniobras comienzan a cobrar mayor importancia, si bien esto es en ocasiones como ya se ha comentado por el funcionamiento del algoritmo (siendo algunas MEP prácticamente nulas o pequeñas correcciones finales). Figura 5.27 Zoom a Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025. Ante la duda de si el Parámetro de Control relacionado con la tolerancia pudiera estar jugando un papel perjudicial en esta simulación por ser del mismo orden de magnitud que el óptimo, se decide repetir el Frente de Pareto bajando la tolerancia (reduciendo el número de puntos calculados, sólo como comprobación) para corroborar que los resultados no están siendo malinterpretados debido a un mal funcionamiento del algoritmo genético. Figura 5.28 Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, con TOL =1.
82 Capítulo 5. Simulaciones y resultados En la Figura 5.28 se observa que la disminución en el Parámetro de Tolerancia no ha afectado a la forma que toma la curva definida por las soluciones no dominadas, a excepción del último tramo donde se han encontrado algunas soluciones mejores 9 . Los tiempos sin embargo se han casi duplicado, de unas 6 horas de procesamiento para la Figura 5.27 a 10 horas en generar el segundo con menor tolerancia (habiéndose evaluado 50 casos menos). No ha de extrañar la diferencia en el número de MEP para algunos tiempos entre ambos Frentes de Pareto, los tiempos de cálculo no son exactamente los mismos, y aún si lo fueran se recuerda que la solución no es fija para un algoritmo metaheurístico (además de la posibilidad de una solución con más maniobras tener alguna de ellas de valor nulo). Demostrada que la forma del Frente de Pareto no es debida a un desajuste en el algoritmo, se procede a analizar algunas soluciones, que se tomarán de la Figura 5.27. En ambos frentes se observa tras una tendencia creciente con 1 MEP un segundo mínimo en la curva, correspondiente con TV=846 días (y ∆VT=5.8931 km/s) y en la Figura 5.27 para 2 MEP. La solución se corresponde con una transferencia directa hasta la órbita planetaria de Marte, añadiendo un phasing para alcanzar luego en este mismo punto al planeta destino. Se puede apreciar que la maniobra de phasing no es excesivamente costosa en este caso y reduce el impulso necesario de establecimiento en la órbita de aparcamiento final si se compara con la Figura 5.26, con lo que el incremento total de impulso necesario respecto a esta no es excesivo. Figura 5.29 Solución con 2 MEP para TV=846 JD. Salida el 29 de Julio de 2024. Otra solución interesante para analizar es la que se corresponde con lo que aparenta ser un máximo local del Frente de Pareto. Se trata de la solución para TV=556 días con 1 MEP, que consigue ∆VT=6.7126 km/s. Se corresponde con una transferencia bielíptica, que consigue llegar al planeta destino en dos arcos de medias elipses. Este tipo de transferencias es óptima en un problema plano de transferencia entre dos órbitas circulares para una relación entre sus radios mayor a 15.58 [ 12 ]. En esta misión donde dicha relación no se cumple, puede comprobarse que no supera a la transferencia 9 La tolerancia no es directamente el error que se comete con respecto al óptimo, se recuerda que es un parámetro que interviene en la decisión de parada del algoritmo y está relacionado con los valores que se van alcanzando en las generaciones para la Función Objetivo de todos los individuos, por lo que depende sustancialmente también del número de población.
5.4 Misión Tierra-Marte 83 directa similar a tipo Hohmann, pero que sí resulta ser la mejor opción para una configuración tan perjudicial como la que se observa en la Figura 5.30, estando los planetas en la salida y en la llegada prácticamente alineados (y no opuestos respecto al Sol), siendo para este tiempo de vuelo aproximadamente donde se alcanzan valores máximos en ∆VT para la misión sin MEP (Figura 5.24). Figura 5.30 Solución con 1 MEP para TV=556 JD. Salida el 27 de Agosto de 2024. Para terminar, se muestra la solución correspondiente a TV=1012 días con 4 MEP, interesante por mostrar la utilidad de un gran número de maniobras en este caso ninguna de ellas nulas pero todas de pequeño valor, logrando un resultado final de ∆VT=6.2078 km/s. Figura 5.31 Solución con 4 MEP para TV=1012 JD. Salida el 21 de Agosto de 2024.
84 Capítulo 5. Simulaciones y resultados Tras los análisis se puede concluir en que la cercanía de ambos planetas hace de las soluciones sin impulsos intermedios una opción muy restringida, delimitada a unos tiempos de vuelo y fechas muy concretos. Las MEP mejoran notablemente el gasto de combustible en el resto de situaciones aunque no consigue mantener el mínimo en el mismo entorno al igual que ocurría en las dos simulaciones anteriores. Más detalles sobre esta misión se muestran para el caso tridimensional en la Sección 5.5. 5.4.1 Análisis de explotación en generación de Frente de Pareto (Tierra-Marte) Para cerrar la sección, se representa de nuevo una gráfica comparativa entre la solución encontrada por el algoritmo y la refinada por @fmincon . En esta sección se aprovecha además que se ha realizado el Frente de Pareto con dos valores del Parámetro de Control TOL para representar ambos casos y comparar y observar una vez más el efecto de este parámetro. En línea continua pueden comprobarse los valores obtenidos para el ajuste de Parámetros de Control con el que se ha trabajado en las simulaciones anteriores (4.2). Las diferencias para este caso entre las soluciones alcanzadas por el algoritmo y las refinadas se mantienen en los mismos márgenes que para las simulaciones anteriores, validando el ajuste aún para esta misión donde el entorno del óptimo de ∆VTes menor. En línea discontinua pueden verse las obtenidas tras bajar TOL a 1. Como puede comprobarse la explotación mejora notablemente (a costa de un alto incremento de tiempo de procesamiento), disminuyendo en todos los casos las diferencias entre las soluciones encontradas por el algoritmo y las refinadas a posteriori. También recalcar que en esta simulación el refinador no ha empeorado ninguna de las soluciones encontradas por el algoritmo. Figura 5.32 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Marte). Aunque este tipo de sucesos pueda parecer una desventaja del algoritmo, la necesidad de realizar comprobaciones cambiando algún parámetro, no es sino una ventaja: la flexibilidad del método. El óptimo de un problema es desconocido antes de su resolución (lógicamente), por lo que una característica beneficiosa de los algoritmos genéticos es la posibilidad de buscarlo bajo distintas opciones de ajuste como se ha realizado en esta sección para corroborar los resultados.
5.5 Ampliación a caso tridimensional 85 5.5 Ampliación a caso tridimensional Se procede ahora a plantear y resolver una misión considerando las órbitas reales de los planetas, es decir, en un caso tridimensional. Se ha decidido escoger para este estudio una misión Tierra-Marte, que ilustra adecuadamente las diferencias con el problema plano y es además la más atractiva a analizar por los planes de futuro cercano relacionados con este planeta. La configuración de misión planteada es similar a la del caso plano: • Fechas de lanzamientos comprendidas entre el 1 de Enero de 2024 y 1 de Enero de 2025. El objetivo es encontrar soluciones viables en cualquier época del año, por lo que el rango es adecuado. • Órbita de aparcamiento inicial terrestre a una altura de 100 km. Se supone con la inclinación adecuada para que la órbita de escape determinada en la solución óptima sea coplanaria a esta. • Órbita final a un radio de 6779 km (dos veces el radio de Marte). Al igual que con la órbita inicial, se asume que la inclinación de esta órbita final es la misma que la de la hipérbola de llegada al planeta. • Tiempos de Vuelo para la generación del Frente de Pareto que comprenden desde un 100 días hasta los 3 años. Con la misma configuración establecida que en el caso plano es posible apreciar los cambios que una inclinación media de diferencia entre ambas órbitas planetarias de tan solo unos 1.85º es capaz de producir. Antes de resolver el problema y generar el Frente de Pareto, conviene comprobar la funcionalidad de los programas para tres dimensiones de forma similar a como se realizó para el caso bidimensional. Comenzando por la representación del Porkchop plot para las fechas de lanzamiento y tiempos de vuelo establecidos: Figura 5.33 Porkchop plot de misión Tierra-Marte, representando C3T. En la Figura 5.33 puede comprobarse que, si bien la disposición de las zonas de coste asumible con una maniobra se encuentran en torno a las mismas fechas y con formas similares al caso bidimensional (Figura 5.25), aparecen dos mínimos (el valor especificado en la Figura 5.33 es orientativo) relacionados con las transferencias Tipo I yTipo II debido a la diferencia de inclinación de las órbitas, como puede ocurrir en casos reales (por ejemplo en la Figura 5.1). Es interesante ver
86 Capítulo 5. Simulaciones y resultados qué caracteriza estos dos tipos de transferencia, por ello se procede a su búsqueda con el algoritmo genético, comprobando así su correcto funcionamiento junto al del resto de programas10. Figura 5.34 Transferencia Tipo I, con fecha de salida el 29 de Octubre de 2024, TV=270 días. Figura 5.35 Transferencia Tipo II, con fecha de salida el 29 de Septiembre de 2024, TV=350 días. Como se puede comprobar, la proyección en el plano es casi idéntica, sin embargo no lo son sus costes ni la inclinación del plano de la órbita de transferencia. Se puede apreciar una diferencia de entorno a 1 km/s entre ambas soluciones. Este suceso es natural, y dependiendo de las fechas de estudio un tipo supera en rentabilidad al otro, siendo en el año de estudio más económico el Tipo II (solución lenta). Es importante aclarar que se está suponiendo en las representaciones anteriores órbita inicial y de aparcamiento final acorde a la configuración de misión, es decir, coplanarias a la hipérbola de salida y de llegada respectivamente, lo que hace que el mínimo en ∆VT se produzca 10 En las representaciones de soluciones la vista tridimensional no mantiene la escala para que se aprecien mejor las diferencias de inclinación entre las distintas órbitas. Su objetivo es ilustrar diferencias.
5.5 Ampliación a caso tridimensional 87 donde se producen los de C3T . De no ser así, la inclinación de llegada jugaría un papel fundamental en la determinación del óptimo. Verificado el funcionamiento de los programas que intervienen, se procede a obtener el Frente de Pareto para la misión. Aún reduciendo a 20 Tiempos Totales de Vuelo evaluados, el tiempo de procesado asciende a las 17 horas y media. Como ya se ha comentado en alguna ocasión, liberar la hipótesis bidimensional no sólo aumenta el número de variables sino también la complejidad del problema, haciendo que el cambio en cualquiera de estas variables tenga mayor efecto sobre la Función Objetivo. Figura 5.36 Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, problema tridimensional. De nuevo y también para el caso tridimensional, se observa una tendencia para las misiones sin impulsos intermedios, que suponen las soluciones óptimas para tiempos cortos. A partir de los 450 días aproximadamente este tipo de soluciones está dominada por aquellas que contienen Maniobras de Espacio Profundo. Figura 5.37 Zoom a Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, problema tridimensional.
88 Capítulo 5. Simulaciones y resultados La forma de la curva de soluciones no dominadas obtenida se asemeja en cierto modo la que aparece en el problema plano, con dos diferencias fundamentales: el coste en general aumenta para toda la curva, y la zona de soluciones MEP óptimas se mantiene más plana en torno a un mismo valor de ∆VT . La primera de las diferencias tiene rápida explicación, reside en la diferencia de planos orbitales entre los dos planetas involucrados. La necesidad de utilizar maniobras para cambiar de plano no solo encarece la misión sino hace menos aprovechable la velocidad del planeta de salida en algunos casos, como los que no presentan impulsos intermedios. Como hasta ahora, se procede a analizar algunas de las soluciones no dominadas para determinar qué hace que el valor de ∆VT oscile entorno a un valor de unos 8-9 km/s para tiempos superiores a los 500 días. Comenzando por analizar la solución con 1 MEP para TV=519 días, en la Figura 5.38 puede apreciarse una solución interesante. Figura 5.38 Solución con 1 MEP para TV=519 JD. Salida el 2 de Marzo de 2024. Si se compara la solución con la de las Figuras 5.34 y 5.35, puede comprobarse que el impulso intermedio trata de provocar una situación similar a una transferencia Tipo II. El primer impulso (el que se realiza en la órbita de aparcamiento inicial) tiene como único objetivo escapar de la esfera de influencia terrestre mediante una parábola, permaneciendo después en la órbita del planeta hasta el momento oportuno de realización de la Maniobra de Espacio Profundo. Se puede demostrar de hecho, que con las condiciones de la órbita inicial, el impulso necesario para escapar de la influencia terrestre mediante una parábola es: ∆V1=r2µ⊕ 6468.14 −rµ⊕ 6468.14 =3.2491 km/s.(5.1) Este primer impulso podría perfectamente haberse ahorrado y combinado con la Maniobra de Espacio Profundo dando lugar a una situación similar a la de la Figura 5.35, pero entonces no se correspondería con una misión del tiempo de vuelo designado. Este resultado es un indicativo de la ventaja de esperar a lanzar en la fecha oportuna, consiguiendo llegar en la misma fecha final (obsérvese la posición de Marte en ambas soluciones) a un coste menor. Esta solución se repite en varios casos de las soluciones no dominadas, en ocasiones cuando no es posible cuadrar bien el primer impulso intermedio (la solución también aparece para más de una MEP no nula) una ligera modificación como la del caso que se aprecia en la Figura 5.39 resulta más efectiva.
5.5 Ampliación a caso tridimensional 89 Figura 5.39 Solución con 1 MEP para TV=624 JD. Salida el 1 de Enero de 2024. No todas las soluciones no dominadas se corresponden sin embargo con esta forma de trayectoria, no posible cumplir con facilidad para tiempos superiores. Para TV=728 días se observa en la solución no dominada de 3 MEP un cambio de plano que se ejecuta en el segundo impulso intermedio, más alejado del Sol para abaratar dicha maniobra (Figura 5.40). Figura 5.40 Solución con 3 MEP para TV=728 JD. Salida el 8 de Marzo de 2024. Para tiempos aún más elevados esta solución se mezcla con algún phasing de espera en alguna de las órbitas de transferencia. Es el caso de la solución de 2 MEP para TV=938 días, que presenta el mayor valor de las soluciones no dominadas con ∆VT=9.5309 km/s (Figura 5.41); o la solución para TV=1043 días con 4 MEP (Figura 5.42).
96 Capítulo 5. Simulaciones y resultados La cantidad de datos se debe a que en este artículo el procedimiento de optimización abordado es distinto, y requiere de un mallado en los tiempos de vuelo y barrido sobre el mismo. Lanzando el algoritmo genético creado en conjunto al refinador, se obtiene el Frente de Pareto de la Figura 5.52. Figura 5.52 Optimización Global de una misión Tierra-Marte. Si se comparan las Figuras 5.52 y 5.51 pueden comprobarse las similitudes, pero también diferencias en el tramo final del Frente de Pareto. Estas diferencias no deben sorprender, se deben principalmente a la diferencia entre las órbitas iniciales y finales (principalmente por la no necesidad de realizar un cambio de plano en la llegada), sumadas a la libertad en la simulación realizada de elegir cualquier Tiempo de Vuelo Total en el rango permitido. En [ 6 ], los autores expresan que Because we search in a grid, instead of randomly evaluating point, the points on Figure 5.51 seem actually ordered. There is however little chance we miss good points in this search space according to the smooth shape of the Pareto front. (Dado que hemos buscado en un espacio discreto, en lugar de evaluar aleatoriamente puntos, los datos en la Figura 5.51 parecen verdaderamente ordenados. Hay sin embargo poca probabilidad de haber perdido buenas soluciones en el espacio de búsqueda acorde a la forma suave del Frente de Pareto.). No obstante, durante el desarrollo de este trabajo ya se ha demostrado que el Tiempo de Vuelo juega un papel crucial en el coste total de la misión, aunque pocas diferencias en este tiempo no supongan cambios en el resultado del orden del óptimo, sí que pueden provocarlo del orden de las diferencias apreciadas. Pasando a analizar qué soluciones componen este nuevo Frente de Pareto, en la Figura 5.53 puede apreciarse que durante un rango considerable de fechas las soluciones sin MEP dominan, pero que durante las primeras fechas quedan muy por encima de las soluciones de 2 MEP. Durante este tramo también las soluciones de un único impulso intermedio mejoran a las transferencias sin MEP.
5.6 Influencia de la fecha de salida. Comparación con otros autores 97 Figura 5.53 Frente de Pareto enfrentando Fecha de Salida y Coste Total (Tierra-Marte, 3D). Para dar finalizar el capítulo de simulaciones, resulta interesante representar el mismo Frente de Pareto enfrentando de nuevo Fecha de Lanzamiento y Tiempo Total de Vuelo. Figura 5.54 Frente de Pareto enfrentando Fecha de Salida y Tiempo de Vuelo Total (Tierra-Marte, 3D). En la Figura 5.54 se puede comprobar que las opciones que minimizan costes con el uso de MEP alargan los tiempos de misión, resultado que concuerda con las simulaciones con Fecha de Salida libre y Tiempo de Vuelo fijado.
6 Conclusiones y trabajo futuro 6.1 Conclusiones Tras el desarrollo del trabajo cabe reflexionar sobre los dos intereses principales del mismo: la utilidad de los algoritmos genéticos en el ámbito aplicado, y la influencia de las Maniobras de Espacio Profundo en las trayectorias interplanetarias. Respecto al primer punto, no cabe duda que el uso de algoritmos genéticos ofrece buenos resultados, que sumado a un método determinista auxiliar para mejorar la explotación supera o como mínimo iguala a los métodos de optimización tradicionales. En cuanto a la cuestión de si merece la pena en cuanto a tiempos de procesamiento y complejidad de programación la respuesta puede conllevar controversia. Si bien el método parece complejo en su descripción, su programación en un lenguaje de alto nivel como puede ser MATLAB es asequible partiendo de unos conocimientos básicos, y su gran adaptabilidad y robustez rentabilizan el esfuerzo. Los tiempos de cómputo, por otro lado, es cierto que son elevados, y otros métodos pueden llegar a ofrecer los mismos resultados más rápidamente. La ventaja principal del uso del algoritmo genético reside sin embargo en la no necesidad de conocer a fondo el problema para resolverlo. Mientras la mayoría de estos métodos alternativos requieren de una estimación inicial suficientemente buena que requiere la resolución (al menos parcial) del problema o uno similar, el uso de algoritmos genéticos no requiere más que la Función Objetivo y las restricciones para de forma autónoma alcanzar el óptimo global (como siempre, con la fiabilidad que permite asegurar que se ha encontrado). La herramienta en cuestión que se ha puesto a prueba en este trabajo tiene aún mucho potencial que mostrar, el tiempo es el coste de su uso pero el alcance de la misma es mucho menos limitado que el de los métodos deterministas, permitiendo añadir gran complejidad al problema planteado sin requerir grandes cambios en su manejo (más allá del reajuste de los Parámetros de Control). En cuanto al uso de Maniobras de Espacio Profundo, queda demostrado que su pueden resultar una gran ventaja, pero son situacionales y en ningún caso de los estudiados suponen una mejora en tiempos y costes a la vez a la transferencia sin impulsos intermedios (siempre y cuando, cómo se ha hecho énfasis, la Fecha de Salida no sea una restricción limitante). Suponen sin embargo una solución idónea para la reducción de costes para Tiempos de Vuelo Totales largos, para una gran horquilla de Fechas de Lanzamiento, y para generar trayectorias que pudieran sacar una mayor rentabilidad científica a la misión. Existen otra serie de escenarios capaces de sacar a relucir la ventaja de uso de las MEP que no se han analizado en este trabajo, como por ejemplo el Rendezvous en una misma órbita, una transferencia esquivando obstáculos durante el trayecto, para órbitas de aparcamiento inicial y final con una inclinación determinada, para enlazar planetas y buscar una 99
100 Capítulo 6. Conclusiones y trabajo futuro trayectoria apoyada en el uso de Maniobras Asistidas por Gravedad... En resumen, los objetivos de investigación del proyecto se han cumplido satisfactoriamente, aunque siendo crítico se podría haber esperado una mayor mejora del Problema de Optimización con el uso de MEP y la utilización del algoritmo genético, lo cierto es que los resultados mostrados confirman el buen funcionamiento de los programas creados y planteamiento del problema, y además de ofrecer un abanico de soluciones sobre el que decidir, abren un mundo de posibilidades en la optimización de las trayectorias interplanetarias en esta primera etapa de optimización bajo hipótesis. 6.2 Trabajo futuro Lo que en este trabajo se ha realizado sienta unas bases sólidas sobre las que seguir construyendo en pos de alcanzar una trayectoria (o un mapa de posibilidades de ellas) que cumpla todos los requisitos que se establezcan para ella. Quedan muchas posibilidades a añadir y tener en consideración de cara a lograr este objetivo, en esta línea las propuestas futuras de cara a mejorar el propio Problema de Optimización serían: • Inclusión en el problema de Maniobras Asistidas por Gravedad. Este tipo de maniobras en el que no se ha entrado en detalle durante la memoria (pues quedan fuera del alcance del proyecto) ofrecen posibilidades de mejora en ∆VT mucho mayores a las alcanzadas con Maniobras de Espacio Profundo. Su uso, sin embargo, es limitado y requiere de configuraciones iniciales idóneas, y es en este escenario donde las Maniobras de Espacio Profundo cobran un papel protagonista ampliando la ventana de lanzamiento y las posibilidades existentes. • Mayores restricciones en las órbitas de aparcamiento final e inicial. La propuesta de una misión real mejor definida que las simulaciones que se han hecho para investigar los objetivos del proyecto podrían arrojar nuevos e interesantes resultados, por ejemplo fijando una inclinación concreta en alguna de estas órbitas. • Modelización de la fase de lanzamiento. Es cierto y como ya se ha comentado que este problema constituye otro Problema de Optimización en sí mismo, pero es innegable la necesidad de reducción de costes en esta etapa. Sería interesante poder optimizar para cada situación ambos problemas que están relacionados al unísono. • Consideración del problema con bajo empuje. Los motores eléctricos son ya en la actualidad una alternativa más que prometedora, pero su uso no permite una modelización del problema con impulsos instantáneos. Por otro lado, el algoritmo genético creado si bien ha resultado ser eficiente, robusto, y arrojar buenos resultados, aún tiene margen de mejora. Propuestas futuras para mejorar el programa serían: • Reducción del número de Parámetros de Control. El ajuste puede resultar un proceso largo por la interrelación entre las variables. Siguiendo el camino de NSGA-II , sería interesante reducir el número de parámetros que intervienen en el algoritmo de cara a hacerlo más manejable y permitir un ajuste más rápido ante un cambio de Función Objetivo. • Reducción de tiempos. Si bien se ha logrado un algoritmo rápido en comparación al genético de referencia ( @ga ), aún es posible reducir más los tiempos. El uso de la herramienta parfor de MATLAB puede ser interesante de implementar, pero requeriría cambios en el algoritmo. • Cambiar el algoritmo de resolución del Problema de Lambert utilizado. Se ha comprobado que gran parte del tiempo de procesamiento se consume en la resolución del Problema de Lambert, por lo que escoger otro más rápido es una opción más que conveniente.
Apéndice A Código Algoritmo Genético Desarrollo del algoritmo genético en MALTAB 1%Algoritmo Genetico 2%-------------------------------------------------------------------------- 3function [xsol,Jsol,Summary]=GA_nDSM_TFG(PC,nPop,nvars,Fobj,A,b,UB,LB) 4%PC recoge los parametros de control del algoritmo, son sus opciones 5%-PC(1)=Numero maximo de generaciones 6%-PC(2)=Tolerancia en la solucion optima 7%-PC(3)=E número de elites 8%-PC(4)=Probabilidad de convertirse en progenitor 9%-PC(5)=Probabilidad de mutar 10 %-PC(6)=Parametro eta_c 11 %-PC(7)=t, participantes del torneo 12 %-PC(8)=sigma, varianza de la Normal en Mutacion 13 %-PC(9)=Probabilidad de que mute un gen 14 %-PC(10)=Numero de cruces aleatorios 15 %-------------------------------------------------------------------------- 16 %Declaracion de variables para acelerar proceso: 17 Elite(PC(3),1:nvars+1)=0; 18 19 %Generacion y evaluacion de poblacion aleatoria inicial: 20 %(la ultima columna almacena el valor de la funcion objetivo para el 21 %individuo que conforma la fila) 22 PopX=GenAleatorio(nPop,nvars,A,b,UB,LB,Fobj); 23 Summary=PopX; 24 %Elite inicial: 25 PopN=PopX; 26 for ii=1:PC(3) 27 [~,ind]=min(PopN(:,end)); 28 Elite(ii,:)=PopN(ind,:); 29 PopN(ind,:)=[]; 30 end 31 %-------------------------------------------------------------------------- 32 %Comienza proceso iterativo: 33 N=1; TOL=100; 34 while N<PC(1) && TOL>PC(2) 35 %Cruce: 36 prf=0; %Bandera auxiliar 37 dd=0; %Numero de descendientes 38 PopD(1:nPop,1:nvars+1)=0; %Espacio para descendientes (por rapidez) 39 for ii=1:nPop 40 if rand<=PC(4) 41 if prf==0 42 pro1=PopX(ii,1:nvars); 101
102 Apéndice A. Código Algoritmo Genético 43 prf=1; 44 else 45 pro2=PopX(ii,1:nvars); 46 prf=0; 47 PopD(dd+1:dd+2,1:nvars)=Cruce(pro1,pro2,PC(6),A,b,LB,UB); 48 dd=dd+2; 49 end 50 end 51 end 52 53 %Cruce extra ente aleatorios 54 for jj=1:PC(10) 55 ind1=randi(nPop,1); 56 ind2=randi(nPop,1); 57 while ind2==ind1 58 ind2=randi(nPop,1); 59 end 60 pro1=PopX(ind1,1:nvars); 61 pro2=PopX(ind2,1:nvars); 62 PopD(dd+1:dd+2,1:nvars)=Cruce(pro1,pro2,PC(6),A,b,LB,UB); 63 dd=dd+2; 64 end 65 66 %Mutacion: 67 %Espacio para mutantes provinientes de evaluados(por rapidez) 68 PopM(1:nPop-PC(3),1:nvars+1)=0; 69 PopEval=PopN; %Almacena los individuos que no es necesario volver a evaluar 70 mm=0; %auxiliar para eliminar de los evaluados los que mutan (n de mutados) 71 for ii=1:nPop-PC(3) 72 if rand<=PC(5) 73 PopM(mm+1,1:nvars)=Mutacion(PopN(ii,1:nvars),PC(9),PC(8),A,b,LB,UB); 74 PopEval(ii-mm,:)=[]; 75 mm=mm+1; 76 end 77 end 78 for ii=1:dd 79 if rand<=PC(5) 80 PopD(ii,1:nvars)=Mutacion(PopD(ii,1:nvars),PC(9),PC(8),A,b,LB,UB); 81 end 82 end 83 %Mutacion completa de individuo aleatorio 84 Elegido=randi(dd+nPop-PC(3),1,1); 85 if Elegido<=length(PopEval(:,1)) 86 PopM(mm+1,1:nvars)=Mutacion(PopEval(Elegido,1:nvars),1,PC(8),A,b,LB,UB); 87 mm=mm+1; 88 PopEval(Elegido,:)=[]; 89 elseif Elegido>length(PopEval(:,1)) && Elegido<=nPop-PC(3) 90 PopM(Elegido-length(PopEval(:,1)),1:nvars)=... 91 Mutacion(PopM(Elegido-length(PopEval(:,1)),1:nvars),1,PC(8),A,b,LB,UB); 92 else 93 PopD(Elegido-nPop+PC(3),1:nvars)=... 94 Mutacion(PopD(Elegido-nPop+PC(3),1:nvars),1,PC(8),A,b,LB,UB); 95 end 96 97 %Evaluacion de los nuevos mienbros (todos menos PopEval y Elite): 98 for ii=1:dd 99 PopD(ii,end)=Fobj(PopD(ii,1:nvars)); 100 end 101 for ii=1:mm 102 PopM(ii,end)=Fobj(PopM(ii,1:nvars)); 103 end 104
103 105 %Elitismo: 106 PopT=[PopEval;PopM(1:mm,:);PopD(1:dd,:);Elite]; 107 for ii=1:PC(3) 108 [~,ind]=min(PopT(:,end)); 109 Elite(ii,:)=PopT(ind,:); 110 PopT(ind,:)=[]; 111 end 112 113 %Seleccion: 114 %1-Se realizan nPop-PC(3) torneos sobre PopT 115 for ii=1:nPop-PC(3) 116 [PopN(ii,1:nvars+1),ind]=Torneo(PopT,PC(7)); 117 PopT(ind,:)=[]; 118 end 119 120 %2-Ordenados, favorece cruce entre mejores 121 [~,s]=sort(PopN(:,end));B=PopN(s,:);PopX=[Elite;B]; 122 123 [Jsol,indsol]=min(Elite(:,end)); 124 xsol=Elite(indsol,1:nvars); 125 TOL=sum(PopX(1:nPop-round(nPop/5),end))/(nPop-round(nPop/5))-Jsol; 126 % if nvars>2 127 % figure(1) 128 % hold on; 129 % plot(N,Jsol,’or’,’MarkerSize’,4); 130 % plot(N,sum(PopX(1:nPop-round(nPop/5),end))/(nPop-round(nPop/5)),... 131 % ’ok’,’MarkerSize’,4); 132 % hold off; 133 % Summary=[Summary;[N,zeros(1,nvars-3),dd,0,mm];PopX]; 134 % end 135 N=N+1; 136 %Ayuda a la convergencia en etapa final 137 if abs(TOL-PC(2))<=1.5*PC(2) || N>PC(1)*4/5 138 PC(7)= 7; 139 end 140 end 141 end 142 %-------------------------------------------------------------------------- 143 %-------------------------------------------------------------------------- 144 %GENERACION DE POBLACION ALEATORIA: 145 %-Se asume que todas las variables estan acotadas entre maximos y minimos 146 function PopX=GenAleatorio(nPop,nvars,A,b,UB,LB,Fobj) 147 PopX(1:nPop,1:nvars+1)=0; 148 for ii=1:nPop 149 PopX(ii,1:nvars)=LB+(UB-LB)*diag(rand(nvars,1))’; 150 %Comprobacion dentro de restricciones 151 while any(A*PopX(ii,1:nvars)’>=b)==1 152 PopX(ii,1:nvars)=LB+(UB-LB)*diag(rand(nvars,1))’; 153 end 154 %Evaluacion con funcion objetivo 155 PopX(ii,end)=Fobj(PopX(ii,1:nvars)); 156 end 157 end 158 %-------------------------------------------------------------------------- 159 %-------------------------------------------------------------------------- 160 %CRUCE 161 function Desc = Cruce(pro1,pro2,etac,A,b,LB,UB) 162 Desc=2*UB(1)*ones(2,length(pro1)); %Fuerza entrar a iterar 163 cont=0; 164 % nvars=length(pro1); 165 while (any(A*Desc(1,:)’>=b)==1 || any(A*Desc(2,:)’>=b)==1 || ... 166 any(Desc(1,:)<LB)==1 || any(Desc(2,:)<LB)==1 || any(Desc(1,:)>UB)==1 ||...
104 Apéndice A. Código Algoritmo Genético 167 any(Desc(2,:)>UB)==1) && cont<1000 168 %- 169 u=rand; 170 if u<=0.5 171 beta=(2*u)^(1/(etac+1)); 172 else 173 beta=(1/2/(1-u))^(1/(etac+1)); 174 end 175 Desc=[0.5*((1+beta)*pro1+(1-beta)*pro2);... 176 0.5*((1-beta)*pro1+(1+beta)*pro2)]; 177 cont=cont+1; 178 end 179 if cont==1000 %No se ha encontrado descendiente 180 %Igual a progenitores 181 beta=1; 182 Desc=[0.5*((1+beta)*pro1+(1-beta)*pro2);... 183 0.5*((1-beta)*pro1+(1+beta)*pro2)]; 184 end 185 end 186 %-------------------------------------------------------------------------- 187 %-------------------------------------------------------------------------- 188 %MUTACION 189 function Mutado=Mutacion(Indiv,pm,sigma,A,b,LB,UB) 190 mut=0; %Variable bandera, se obliga a al menos una mutacion 191 Mutado=Indiv; 192 for ii=1:length(Indiv) 193 u=rand; 194 if pm>=u 195 n=normrnd(0,(UB(ii)-LB(ii))/sigma); 196 Mutado(ii)=Indiv(ii)+n; 197 cont=0; 198 while (Mutado(ii)>UB(ii) || Indiv(ii)+n<LB(ii) || ... 199 any(A*Mutado’>=b)==1) && cont<1000 200 % n=normrnd(0,sigma); 201 n=normrnd(0,(UB(ii)-LB(ii))/sigma); 202 cont=cont+1; 203 Mutado(ii)=Indiv(ii)+n; 204 end 205 if cont==1000 206 Mutado(ii)=Indiv(ii); 207 end 208 mut=1; 209 end 210 end 211 if mut==0 212 cont=0; 213 ii=randi(length(Indiv),1,1); 214 n=normrnd(0,(UB(ii)-LB(ii))/sigma); 215 Mutado(ii)=Indiv(ii)+n; 216 while (Mutado(ii)>UB(ii) || Indiv(ii)+n<LB(ii) || ... 217 any(A*Mutado’>=b)==1) && cont<1000 218 n=normrnd(0,(UB(ii)-LB(ii))/sigma); 219 Mutado(ii)=Indiv(ii)+n; 220 cont=cont+1; 221 end 222 if cont==1000 223 Mutado(ii)=Indiv(ii); 224 end 225 end 226 end 227 %-------------------------------------------------------------------------- 228 %--------------------------------------------------------------------------
105 229 %SELECCION POR TORNEO 230 function [Sel,ind]=Torneo(PopT,t) 231 [m,~]=size(PopT); 232 Candidatos=sort(randi(m,1,t)); 233 for ii=1:t-1 234 if Candidatos(ii)==Candidatos(ii+1) 235 Candidatos(ii+1)=Candidatos(ii+1)+1; 236 end 237 end 238 while Candidatos(t)>m || ... 239 any(Candidatos(t)*ones(1,t-1)==Candidatos(1:t-1))==1 240 Candidatos(t)= randi(m,1,1); 241 Candidatos=sort(Candidatos); 242 end 243 [~,ind]=min(PopT(Candidatos,end)); 244 Sel=PopT(ind,:); 245 end
112 Capítulo B. Herramientas de modelización del problema 183 %de salida|impulsos, excepto el impulso n|IncrVtotal 184 185 if ga_opt~=0 186 %Optimos sacados por el algoritmo genetico: 187 Datosga(1:AA*(n+1-n_inf),1:3*n+7)=NaN; 188 for ii=1:(n+1-n_inf)*AA 189 [Datosga(ii,3),Datosga(ii,2),Datosga(ii,1),~]=Jdate_inv(Xga(ii,1)); 190 [Datosga(ii,6),Datosga(ii,5),Datosga(ii,4),~]=Jdate_inv(Xga(ii,1)+Tv); 191 for kk=2:n+1 192 Datosga(ii,kk+5)=sum(Xga(ii,2:kk)); 193 end 194 end 195 Datosga(:,n+7:3*n+6)=Xga(:,n+2:3*n+1); 196 Datosga(:,3*n+7)=Jga’; 197 end 198 199 %Los optimos refinados 200 Datos(1:AA*(n+1-n_inf),1:3*n+7)=NaN; 201 for ii=1:(n+1-n_inf)*AA 202 [Datos(ii,3),Datos(ii,2),Datos(ii,1),~]=Jdate_inv(X(ii,1)); 203 [Datos(ii,6),Datos(ii,5),Datos(ii,4),~]=Jdate_inv(X(ii,1)+Tv); 204 for kk=2:n+1 205 Datos(ii,kk+5)=sum(X(ii,2:kk)); 206 end 207 end 208 Datos(:,n+7:3*n+6)=X(:,n+2:3*n+1); 209 Datos(:,3*n+7)=J’; 210 toc 211 212 %-------------------------------------------------------------------------- 213 %Creacion del Frente de Pareto con las soluciones refinadas: 214 if AA>1 215 figure;hold on; 216 Dibn=[’*’,’o’,’+’,’^’,’x’]; 217 for aa=1:AA 218 plot(TV(aa),Jmin(aa),Dibn(Ind(aa)),’color’,’r’,’MarkerSize’,8); 219 if Ind(aa)==1 220 for jj=2:n-n_inf+1 221 plot(TV(aa),J(jj+(aa-1)*(n-n_inf+1)),Dibn(jj),’color’,’k’,’MarkerSize’,6); 222 end 223 else 224 for jj=[1:Ind(aa)-1,Ind(aa)+1:n-n_inf+1] 225 plot(TV(aa),J(jj+(aa-1)*(n-n_inf+1)),Dibn(jj),’color’,’k’,’MarkerSize’,6); 226 end 227 end 228 end 229 end 230 %-------------------------------------------------------------------------- 231 %-------------------------------------------------------------------------- 232 233 %-------------------------------------------------------------------------- 234 %Funcion Objetivo 235 236 function J = Fobjetivo_nDSM(x) 237 global hi; global planet; global Rff; global Tv; 238 n=(length(x)-1)/3; 239 mu=[22032 324859 398600.4 42828 126711995.4 37931187 5793939 ... 240 6836529 871]; 241 %x(1) es fecha de salida 242 %x(2:n+1) son las fracciones de tiempo de vuelo entre maniobras 243 %x(n+2:3*n+1) son las componentes de velocidad inicial y de las n-1 244 %maniobras
B.3 Código del Problema de Optimización y Función Objetivo 113 245 IV(1:n+2)=NaN;%vector de incrementos de velocidad 246 %Situacion inicial y final: 247 [r0,v0] = efemeride(3,x(1)); 248 [rf,vf] = efemeride(planet,x(1)+Tv); 249 %-------------------------------------------------------------------------- 250 %Caso sin maniobras: 251 %-------------------------------------------------------------------------- 252 if n==0 253 [v1,v2,flag]=LambertV2(Tv,r0,rf,10000); 254 [~,~,~,~,e,~,~,~,~]=Elementos(r0,v1,1); 255 if imag(v1(1))~=0 || imag(v1(2))~=0 || imag(v2(1))~=0 || imag(v2(2))~=0 256 IncrVt=NaN; 257 elseif isnan(e)==1 || flag==10000 258 IncrVt=NaN; 259 else 260 C3i=norm(v1-v0)^2; 261 IV(1)=sqrt(2*398600.4/(6378.14+hi)+29.7847^2*C3i)-sqrt(398600.4/(6378.14+hi)); 262 C3f=norm(v2-vf)^2; 263 IV(2)=sqrt(2*mu(planet)/Rff+29.7847^2*C3f)-sqrt(mu(planet)/Rff); 264 IncrVt=abs(IV(1))+abs(IV(2)); 265 end 266 %-------------------------------------------------------------------------- 267 %Caso con maniobras: 268 %-------------------------------------------------------------------------- 269 else 270 %Actualmente interpretacion de v1 como velocidad de escape de la tierra 271 gamma=x(n+3); 272 C3i=x(n+2)^2; 273 v1=x(n+2)*v0./norm(v0); %Esta seria la velocidad tangente 274 C=[cos(gamma),sin(gamma),0;-sin(gamma),cos(gamma),0;0,0,1]; 275 v1=(C*v1’)’; 276 vi=v1+v0; %SIEMPRE EN UNIDADES CANONICAS 277 IV(1)=sqrt(2*398600.4/(6378.14+hi)+29.7847^2*C3i)-sqrt(398600.4/(6378.14+hi)); 278 %A partir de este punto ri se actualiza como la posicion al inicio de conica 279 %y vi como su velocidad en este punto. 280 ri=r0; 281 for kk=1:n 282 [RAAN,w,i,p1,e1,theta1,w_,u,lambda1]=Elementos(ri,vi,1); 283 %---------------------------------------------------- 284 %Calculo de la posicion tras el tiempo de vuelo hasta DSM 285 if theta1>=0 286 T1=tvuelo(0,theta1,p1,e1); 287 T=x(1+kk)*Tv+T1; 288 theta2=LeyHoraria(p1,e1,T); 289 elseif theta1<0 290 T1=tvuelo(theta1,0,p1,e1); 291 if T1<x(1+kk)*Tv 292 T=x(1+kk)*Tv-T1; 293 theta2=LeyHoraria(p1,e1,T); 294 else 295 T=T1-x(1+kk)*Tv; %en este caso el angulo obtenido sera -theta2 296 theta2=LeyHoraria(p1,e1,T); 297 theta2=2*pi-theta2; 298 end 299 end 300 lambda2=[]; 301 if isempty(theta1)==1 %Es circular 302 Incrlambda=LeyHoraria(p1,e1,x(1+kk)*Tv); 303 lambda2=lambda1+Incrlambda; 304 lambda2=lambda2-2*pi*floor(lambda2/2/pi); %por si diera mas de una vuelta 305 theta2=[]; 306 end
114 Capítulo B. Herramientas de modelización del problema 307 [ri,v2]=posvel(RAAN,w,i,p1,e1,theta2,w_,u,lambda2,1); 308 %si aun no se ha llegado a la ultima maniobra: 309 if kk<n 310 vi=v2+[x(n+2+2*kk),x(n+3+2*kk),0]./29.7847; %se le suma la DSM en UVsol 311 IV(kk+1)=norm([x(n+2+2*kk),x(n+3+2*kk),0]); %en km/s 312 end 313 end 314 %La ultima DSM esta fijada por el problema de Lambert para asegurar llegada 315 %El tiempo de vuelo de esta maniobra es el restante hasta la fecha de 316 %llegada, se pone en absoluto para que no de fallos el programa, en las 317 %constraint se comprueba que el tiempo hasta entonces no ha superado el 318 %total 319 Tfin=abs(1-sum(x(2:n+1)))*Tv; 320 [v3,v4,flag]=LambertV2(Tfin,ri,rf,5000); 321 if flag<5000 322 IV(n+1)=29.7847*norm(v3-v2); %incremento de velocidad de la DSM, en km/s 323 else 324 IV(n+1)=Inf;%no ha encontrado solucion 325 end 326 %Por ultimo la llegada a la orbita de aparcamiento 327 C3f=norm(v4-vf)^2; %en UVsol 328 IV(n+2)=abs(sqrt(2*mu(planet)/Rff+29.7847^2*C3f)-sqrt(mu(planet)/Rff)); 329 IncrVt=sum(IV); 330 %-------------------------------------------------------------------------- 331 clc 332 end 333 J = IncrVt; %en km/s 334 end 335 %-------------------------------------------------------------------------- 336 %--------------------------------------------------------------------------
Índice de Figuras 1.1 Misiones de espacio profundo de mayor relevancia dentro del Sistema Solar, anteriores al año 2003 1 1.2 Misiones de espacio profundo de mayor relevancia dentro del Sistema Solar, posteriores al año 2003 y proyectadas hasta 2018 2 2.1 Sistema de referencia genérico para órbita 8 2.2 Esquema de órbita elíptica en plano perifocal, extraída de [12] 9 2.3 Esquema de órbita parabólica en plano perifocal, extraída de [12] 10 2.4 Esquema de órbita hiperbólica en plano perifocal, extraída de [12] 10 2.5 Definición gráfica de la Anomalía Excéntrica (E) 14 2.6 Representación geométrica del Ángulo Hiperbólico (H) 15 2.7 Propagación en órbita elíptica 15 2.8 Esquema del problema de Lambert, extraída de [2] 17 2.9 Sistema de Referencia Heliocéntrico Eclíptico 19 2.10 Esquema de Maniobra de Espacio Profundo 21 2.11 Velocidad de salida del planeta de salida 22 3.1 Clasificación de problemas según su complejidad 28 3.2 Simulación de tiempos de computación de la función objetivo dada la solución, hasta 50 Maniobras de Espacio Profundo 29 3.3 Diferencia entre exploración y explotación 31 3.4 Diagrama de proceso de Selección por Torneo 33 3.5 Diagrama de proceso de Cruce 35 3.6 Diagrama de proceso de Mutación 36 3.7 Diagrama de flujo del proceso del algoritmo genético 38 3.8 Ejemplo de Pareto Dominancia, extraído de [9] 40 3.9 Ejemplo de Frente de Pareto, extraído de [9] 40 3.10 Ejemplo de distintas formas de Frente de Pareto, extraído de [9] 41 3.11 Ejemplo de Frente de Pareto obtenido por NSGA-II frente a método alternativo (PAES), sobre un problema típico usado para probar optimizadores conocido como ZTD4, extraído de [3] 42 3.12 Esquema de procedimiento de NSGA-II, extraído de [3] 42 3.13 Ejemplo de Frente de Pareto obtenido por NSGA-II ajustado sobre un problema típico usado para probar optimizadores conocido como ZTD4, extraído de [3] 43 4.1 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=50 50 4.2 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=100 51 115
116 Índice de Figuras 4.3 Zoom a mejora del óptimo en procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y X=50 51 4.4 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y N=100. Compárese con la Figura 4.1 53 4.5 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1) y TOL =1. Compárese con la Figura 4.1 54 4.6 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), t=5yE=7. Compárese con la Figura 4.1 55 4.7 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), pc=0.4, ηc=3, y nc=2. Compárese con la Figura 4.1 55 4.8 Procedimiento del algoritmo con los Parámetros de Control iniciales (4.1), pmut =0.15, pm=0.3, y σc=30. Compárese con la Figura 4.1 56 4.9 Procedimiento del algoritmo con los Parámetros de Control finales 57 4.10 Procedimiento del algoritmo en análisis de detalle 58 4.11 Zoom en procedimiento del algoritmo en análisis de detalle 59 4.12 Solución alcanzada por algoritmo genético 60 4.13 Solución alcanzada tras refinado 60 4.14 Procedimiento del algoritmo en verificación del problema 3D 61 4.15 Procedimiento del algoritmo en verificación del problema 3D. Mejor individuo de cada genración 62 5.1 Ejemplo (extraído) de PorkChop plot para misión Tierra-Marte, caso tridimensional 64 5.2 Porkchop plot de misión Tierra-Saturno, representando C3T64 5.3 Solución minimizando ∆VT65 5.4 Solución minimizando C3T65 5.5 Frente de Pareto para una misión Tierra Saturno 66 5.6 Zoom a Frente de Pareto para una misión Tierra Saturno 67 5.7 Frente de Pareto para una misión Tierra Saturno 68 5.8 Tendencia de misiones sin MEP en Frente de Pareto para una misión Tierra Saturno 68 5.9 Zoom a Frente de Pareto para una misión Tierra Saturno con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 69 5.10 Solución con 2 MEP para TV=6935 JD. Salida el 6 de Noviembre de 2024 69 5.11 Solución con 4 MEP para TV=8760 JD. Salida el 3 de Enero de 2024 70 5.12 Solución con 1 MEP para TV=4015 JD. Salida el 3 de Agosto de 2024 71 5.13 Solución con 2 MEP para TV=5840 JD. Salida el 22 de Septiembre de 2024 71 5.14 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Saturno) 72 5.15 Frente de Pareto para una misión Tierra-Mercurio con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 73 5.16 Zoom a Frente de Pareto para una misión Tierra-Mercurio con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 74 5.17 Solución con 1 MEP para TV=1237 JD. Salida el 9 de Junio de 2024 75 5.18 Solución con 4 MEP para TV=1689 JD. Salida el 23 de Mayo de 2024 75 5.19 Solución con 1 MEP para TV=181 JD. Salida el 24 de Mayo de 2024 76 5.20 Solución con 4 MEP para TV=633 JD. Salida el 22 de Mayo de 2024 76 5.21 Solución con 3 MEP para TV=2141 JD. Salida el 23 de Diciembre de 2024 77 5.22 Frente de Pareto para una misión Tierra-Mercurio con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 78 5.23 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Mercurio) 78
Índice de Figuras 117 5.24 Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025. Tendencia de misiones con 0 MEP destacada 79 5.25 Porkchop plot del C3Tpara una misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 80 5.26 Solución con 0 MEP para TV=307 JD. Salida el 3 de Octubre de 2024 80 5.27 Zoom a Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025 81 5.28 Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, con TOL =181 5.29 Solución con 2 MEP para TV=846 JD. Salida el 29 de Julio de 2024 82 5.30 Solución con 1 MEP para TV=556 JD. Salida el 27 de Agosto de 2024 83 5.31 Solución con 4 MEP para TV=1012 JD. Salida el 21 de Agosto de 2024 83 5.32 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Marte) 84 5.33 Porkchop plot de misión Tierra-Marte, representando C3T85 5.34 Transferencia Tipo I, con fecha de salida el 29 de Octubre de 2024, TV=270 días 86 5.35 Transferencia Tipo II, con fecha de salida el 29 de Septiembre de 2024, TV=350 días 86 5.36 Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, problema tridimensional 87 5.37 Zoom a Frente de Pareto para la misión Tierra-Marte con fecha de lanzamiento entre el 1 de Enero de 2024 y el 1 de Enero de 2025, problema tridimensional 87 5.38 Solución con 1 MEP para TV=519 JD. Salida el 2 de Marzo de 2024 88 5.39 Solución con 1 MEP para TV=624 JD. Salida el 1 de Enero de 2024 89 5.40 Solución con 3 MEP para TV=728 JD. Salida el 8 de Marzo de 2024 89 5.41 Solución con 2 MEP para TV=928 JD. Salida el 16 de Julio de 2024 90 5.42 Solución con 4 MEP para TV=1043 JD. Salida el 5 de Marzo de 2024 90 5.43 Diferencias entre soluciones alcanzadas por uso único de algoritmo genético y las refinadas (Tierra-Marte, 3D) 91 5.44 Comparativa entre curva de soluciones no dominadas (Tierra-Marte, 3D) 91 5.45 Comparativa entre Frentes de Pareto (Tierra-Marte, 3D) 92 5.46 Zoom a comparativa entre Frentes de Pareto (Tierra-Marte, 3D) 92 5.47 Frente de Pareto en tres dimensiones, incluyendo fecha de salida (Tierra-Marte, 3D) 93 5.48 Frente de Pareto enfrentando fecha de salida y tiempo de vuelo (Tierra-Marte, 3D) 93 5.49 Frente de Pareto considerando como objetivo fecha de salida (Tierra-Marte, 3D) 94 5.50 Frente de Pareto considerando como objetivo fecha de salida (Tierra-Marte, 3D) 94 5.51 Optimización Global de una misión Tierra-Marte, extraído de [6] 95 5.52 Optimización Global de una misión Tierra-Marte 96 5.53 Frente de Pareto enfrentando Fecha de Salida y Coste Total (Tierra-Marte, 3D) 97 5.54 Frente de Pareto enfrentando Fecha de Salida y Tiempo de Vuelo Total (Tierra-Marte, 3D) 97
Bibliografía [1] F. Cacciatore and C. Toglia, Optimization of Orbital Trajectories using Genetic Algorithms, DEIMOS Space S.L., Madrid, Spain, and University of Rome ’La Sapienza’, Rome, Italy, 2008. [2] Howard D. Curtis, Orbital Mechanics For Engineering Students, ELSEVIER, 2005. [3] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, A fast and elitist multiobjective genetic algorithm: Nsga-ii, IEEE Transactions on Evolutionary Computation 6 (2002), no. 2, 182–197. [4] Kalyanmoy Deb and Amarendra Kumar, Real-coded Genetic Algorithms with Simulated Binary Crossover: Studies on Multimodal and Multiobjective Problems, Department of Mechanical Engineering, Indian Institute of Technology, Kanpur, UP 208 016, India, 1995. [5] Michale R. Garey and David S. Johnson, Computers and Intractability, A Guide to the Theory of NP-Completeness, W.H. FREEMAN AND COMPANY, San Francisco, 1979. [6] Jean-Paul Marmorat Joris T. Olympio, Global Trajectory Optimisation: Can we prune the solution space when considering Deep Space Maneuvers? Final Report, extraído de http://esa.int/act, 2020. [7] Emily Lakdawalla, Planetary Exploration Timelines: A Look Ahead to 2016, December 2015. [8] E.R. Lancaster and R.C.Blanchard, A unified form of Lambert’s Problem, Goddard Space Flight Center, Greenbelt, Md., 1969. [9] Sean Luke, Essentials of Metaheuristics, second ed., Lulu, 2013, Disponible gratis en http://cs.gmu.edu/∼sean/book/metaheuristics/. [10] Jacob A. Englander Matthew A. Vavrina and Donald H. Ellison, Global Optimization of n-Maneuver, High-Thrust Trajectories using Direct Multiple Shooting, NASA. [11] E.M. Standish, Keplerian elements for approximate positions of the major planets, Solar System Dynamics Group, JPL/Caltech, 2006. [12] Rafael Vázquez Valenzuela, Astronáutica y Vehículos Espaciales, Apuntes de clase (ETSI Sevilla), extraídos de http://aero.us.es/astro/desc.html, 2019-2020. [13] Xin-She Yang, Nature-Inspired Metaheuristic Algorithms, Middlesex University, UK, 2010. 119