Machine Learning aplicado a la propagación de alta precisión de trayectorias de satélites

TRABAJO FIN DE GRADO Machine Learning aplicado a la propagaci´on de alta precisi´on de trayectorias de sat´elites Realizado por Jos´e Calder´on Valdivia Para la obtenci´on del t´ıtulo de Doble Grado en Ingenier´ıa Inform´atica - Tecnolog´ıas Inform´aticas y Matem´aticas Dirigido por Daniel Ayala Hern´andez Luis Valencia Cabrera Realizado en el departamento de Lenguajes y Sistemas Inform´aticos Convocatoria de junio, curso 2022/23 Agradecimientos Quiero agradecer a mi familia por apoyarme en todo momento y haber estado ah´ı siempre que lo he necesitado; sin ellos no habr´ıa sido posible. A mis amigos, cuya compa˜n´ıa ha sido un pilar fundamental durante estos a˜nos de carrera. A mis profesores por todo lo que me han ense˜nado. A todo el grupo de investigaci´on DEAL, que desde el primer d´ıa me acogieron como a uno m´as. A Rafael Ayala, autor del paquete asteRisk, por resolver muchas dudas que surgieron durante el desarrollo del proyecto. Y especialmente quiero agradecer a mis tutores, Daniel Ayala y Luis Valencia, su ayuda a lo largo de todo el proyecto. i Resumen La predicci´on precisa de las posiciones de los sat´elites es fundamental para una amplia gama de aplicaciones, que van desde las comunicaciones y navegaci´on por sat´elite hasta la teledetecci´on y la planificaci´on de misiones espaciales. La capacidad de predecir con precisi´on las trayectorias de los sat´elites proporciona informaci´on cr´ıtica para el funcionamiento seguro y eficiente de estos sistemas espaciales. Sin embargo, debido a las diversas fuerzas y perturbaciones que afectan a los sat´elites en ´orbita, lograr una predicci´on precisa enfrenta desaf´ıos t´ecnicos y cient´ıficos significativos. Para la modelizaci´on de trayectorias existen principalmente dos tipos de modelos: los anal´ıticos, como SGP4 y SDP4; y los de alta precisi´on, como HPOP. Como su nombre indica, los de la segunda familia aportan m´as fiabilidad; sin embargo, el coste computacional es demasiado elevado para realizar predicciones a medio o largo plazo. En el presente trabajo se utiliza el paquete de R asteRisk para la creaci´on de un conjunto de datos de entrenamiento y validaci´on y, posteriormente, desarrollar modelos de redes neuronales para corregir las predicciones realizadas por un modelo anal´ıtico y as´ı obtener predicciones precisas en un tiempo reducido. Previamente a la creaci´on del conjunto de datos se desarrolla un m´etodo de optimizaci´on de los par´ametros del modelo SGP4 y SDP4, de modo que en el conjunto de datos se pudieron incluir tanto las predicciones con los par´ametros base como con los par´ametros optimizados. Por tanto, tambi´en se pudo estudiar c´omo variaba la capacidad de predicci´on del error de las redes neuronales seg´un si los par´ametros orbitales de SGP4 y SDP4 hab´ıan sido previamente optimizados o no. Todo esto es posible gracias la gran cantidad de datos disponibles sobre trayectorias de sat´elites en un formato estandarizado, ya sean archivos RINEX o archivos TLE. Palabras clave: aprendizaje autom´atico, redes neuronales, propagaci´on, trayectoria, sat´elite, modelos anal´ıticos, modelos de alta precisi´on ii Abstract Accurate prediction of satellite positions is crucial for a wide range of applications, ranging from satellite communications and navigation to remote sensing and space mission planning. The ability to accurately predict satellite trajectories provides critical information for the safe and efficient operation of these space systems. However, due to the various forces and perturbations affecting satellites in orbit, achieving accurate prediction faces significant technical and scientific challenges. For trajectory modelling, there are mainly two types of models: analytical models, such as SGP4 and SDP4; and high-precision models, such as HPOP. As the name suggests, the latter family provides greater reliability; however, the computational cost is prohibitively high for making medium or long-term predictions. In this work, the R package asteRisk is used to create a training and validation dataset and, subsequently, develop neural network models to correct the predictions made by an analytical model, thus obtain accurate predictions in a reduced time frame. Prior to the creation of the dataset, a parameter optimization method is developed for the SGP4 and SDP4 models, allowing both baseline and optimized parameter predictions to be included in the dataset. Therefore,the study also investigates how the neural networks’ prediction error varies depending on whether the orbital parameters of SGP4 and SDP4 have been previously optimized or not. All of this is made possible by the wealth of available data on satellite trajectories in a standardized format, whether it be RINEX files or TLE files. Keywords: machine learning, neural networks, propagation, trajectory, satellite, analytical models, high-precision models iii ´ Indice general 1. Introducci´on 1 2. Estudio Previo 3 2.1. Objetivos .................................. 3 2.2. Metodolog´ıa................................. 3 2.3. Planificaci´on................................. 4 2.4. Presupuesto................................. 6 3. Estudio Te´orico 8 3.1. Conceptosprevios.............................. 8 3.1.1. Geometr´ıa de la elipse . . . . . . . . . . . . . . . . . . . . . . . 8 3.1.2. Algunas definiciones . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.3. LeyesdeKepler........................... 10 3.1.4. Ecuaci´on de Kepler . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.5. Constelaciones de sat´elites . . . . . . . . . . . . . . . . . . . . . 12 3.2. Sistemas de referencia de coordenadas . . . . . . . . . . . . . . . . . . . 13 3.3. Sistemasdetiempo............................. 14 3.4. Elementosorbitales............................. 16 3.5. Formatos de intercambio de datos de posiciones . . . . . . . . . . . . . 18 3.5.1. TLE (Two-Line Element set) . . . . . . . . . . . . . . . . . . . 18 3.5.2. RINEX (Receiver INdependent EXchange) . . . . . . . . . . . . 18 3.6. Modelos de perturbaciones simplificadas . . . . . . . . . . . . . . . . . 19 3.7. Modelos de propagaci´on de alta precisi´on . . . . . . . . . . . . . . . . . 24 3.7.1. Formulaci´on de Encke . . . . . . . . . . . . . . . . . . . . . . . 25 3.7.2. Formulaci´on de Cowell . . . . . . . . . . . . . . . . . . . . . . . 26 3.7.3. HPOP................................ 26 3.7.4. Resoluci´on num´erica de ecuaciones diferenciales . . . . . . . . . 29 4. An´alisis del problema 33 4.1. Requisitos de informaci´on . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2. Requisitos funcionales . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.3. Requisitos no funcionales . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4. Estadodelarte ............................... 34 5. Dise˜no de la soluci´on 35 5.1. Ideainicial.................................. 35 5.2. Modificaci´on del dise˜no inicial . . . . . . . . . . . . . . . . . . . . . . . 36 5.3. Creaci´on de los modelos . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.4. Evaluaci´on de resultados . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.5. Ampliaci´on de la experimentaci´on . . . . . . . . . . . . . . . . . . . . . 40 6. Problemas encontrados 41 6.1. Ajuste de par´ametros . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 iv 6.2. Representaci´on de ´orbitas . . . . . . . . . . . . . . . . . . . . . . . . . 41 6.3. Repositorios................................. 42 6.3.1. Estaci´on de la Universidad de Sevilla . . . . . . . . . . . . . . . 42 7. Implementaci´on 43 7.1. Herramientas ................................ 43 7.1.1. R................................... 43 7.1.2. asteRisk y asteRiskData . . . . . . . . . . . . . . . . . . . . . . 43 7.1.3. Dplyr ................................ 44 7.1.4. Plotly ................................ 44 7.1.5. Reticulate.............................. 44 7.1.6. Python................................ 44 7.1.7. BeautifulSoup ........................... 45 7.1.8. TensorFlow ............................. 45 7.1.9. CUDAycuDNN .......................... 45 7.2. Desarrollo .................................. 45 7.2.1. Instalaci´on de asteRisk y asteRiskData . . . . . . . . . . . . . . 45 7.2.2. C´alculo de elementos cl´asicos orbitales . . . . . . . . . . . . . . 45 7.2.3. Calcular el error medio de predicci´on . . . . . . . . . . . . . . . 46 7.2.4. Optimizaci´on de par´ametros orbitales . . . . . . . . . . . . . . . 48 7.2.5. Representaci´on de ´orbitas . . . . . . . . . . . . . . . . . . . . . 52 7.2.6. Creaci´on del primer dataset .................... 53 7.2.7. Obtenci´on de datos RINEX . . . . . . . . . . . . . . . . . . . . 53 7.2.8. An´alisis y procesado de los datos . . . . . . . . . . . . . . . . . 55 7.2.9. Creaci´on del dataset con las predicciones . . . . . . . . . . . . . 58 7.2.10. Adici´on de las predicciones con par´ametros sin optimizar . . . . 62 7.2.11. Instalaci´on de TensorFlow y Keras . . . . . . . . . . . . . . . . 64 7.2.12. Instalaci´on de CUDA y cuDNN . . . . . . . . . . . . . . . . . . 64 7.2.13. Implementaci´on de Robust Scaler . . . . . . . . . . . . . . . . . 64 7.2.14. Entrenamiento de los primeros modelos . . . . . . . . . . . . . . 64 7.2.15. Obtenci´on de los datos de los astros . . . . . . . . . . . . . . . . 69 7.2.16. Adici´on de los datos de los astros al dataset ........... 73 7.2.17. Creaci´on de modelos con datos de los astros . . . . . . . . . . . 73 7.2.18. Evaluaci´on de los modelos . . . . . . . . . . . . . . . . . . . . . 75 7.3. Conclusiones................................. 76 8. Pruebas y experimentos 77 8.1. Error medio de las predicciones SGDP4 . . . . . . . . . . . . . . . . . . 77 8.2. Evaluaci´on de los modelos sin datos de los astros . . . . . . . . . . . . 78 8.2.1. Par´ametros orbitales sin optimizar . . . . . . . . . . . . . . . . 78 8.2.2. Par´ametros orbitales optimizados . . . . . . . . . . . . . . . . . 81 8.3. Evaluaci´on de los modelos con datos de los astros . . . . . . . . . . . . 83 8.3.1. Par´ametros orbitales sin optimizar . . . . . . . . . . . . . . . . 84 8.3.2. Par´ametros orbitales optimizados . . . . . . . . . . . . . . . . . 86 8.4. Conclusiones................................. 88 9. Conclusiones 91 v 10.Bibliograf´ıa 92 vi ´ Indice de figuras 3.1. Esquemadeunaelipse ........................... 9 3.2. Situaci´on de una ´orbita . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3. ZonasUTC ................................. 14 3.4. Elementos angulares de una ´orbita . . . . . . . . . . . . . . . . . . . . 17 3.5. Elementos orbitales alternativos . . . . . . . . . . . . . . . . . . . . . . 17 3.6. Formulaci´on de Encke . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.7. Representaci´on gr´afica del m´etodo de Runge-Kutta de cuarto orden . . 30 5.1. Estructura de las redes neuronales implementadas . . . . . . . . . . . . 39 7.1. Estructura de directorios del repositorio de RINEX versi´on 2 de la Junta deAndaluc´ıa ................................ 53 8.1. Evoluci´on del MAE recibiendo los par´ametros base y la posici´on predicha con los par´ametros orbitales sin optimizar . . . . . . . . . . . . . . . . 78 8.2. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial y la posici´on predicha con los par´ametros orbitales sin optimizar . . . . 79 8.3. Evoluci´on del MAE recibiendo los par´ametros base y el movimiento predicho con los par´ametros orbitales sin optimizar . . . . . . . . . . . . . 80 8.4. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial y el movimiento predicho con los par´ametros orbitales sin optimizar . . 80 8.5. Evoluci´on del MAE recibiendo los par´ametros base y la posici´on predicha con los par´ametros orbitales optimizados . . . . . . . . . . . . . . . . . 81 8.6. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial y la posici´on predicha con los par´ametros orbitales optimizados . . . . . 82 8.7. Evoluci´on del MAE recibiendo los par´ametros base y el movimiento predicho con los par´ametros orbitales optimizados . . . . . . . . . . . . . . 82 8.8. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial y el movimiento predicho con los par´ametros orbitales optimizados . . . 83 8.9. Evoluci´on del MAE recibiendo los par´ametros base, los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitalessinoptimizar ........................... 84 8.10. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitales sin optimizar . . . . . . . . . . . . . . . . 85 8.11. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho con los par´ametros orbitales sin optimizar . . . . . . . . . . . 85 8.12. Evoluci´on del MAE recibiendo los par´ametros base, los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitalesoptimizados............................ 86 8.13. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitales optimizados . . . . . . . . . . . . . . . . . 87 vii 8.14. Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho con los par´ametros orbitales optimizados . . . . . . . . . . . . 88 viii Tarea Dedicaci´on estimada Dedicaci´on real Investigaci´on 116 horas 114 horas Recopilaci´on de art´ıculos 2 horas 2 horas Lectura en profundidad de art´ıculos 4 horas 6 horas Fundamentos de las ´orbitas 30 horas 25 horas Modelos de propagaci´on 60 horas 65 horas Resoluci´on num´erica de ec. dif. 20 horas 16 horas Implementaci´on de funciones 45 horas 103 horas Obtenci´on de elementos orbitales 7 horas 5 horas C´alculo del error de predicci´on 3 horas 2 horas Optimizaci´on de par´ametros orbitales 15 horas 60 horas Representaci´on de ´orbitas 3 horas 3 horas Generaci´on de predicciones 10 horas 25 horas Unificaci´on de los datos 5 horas 5 horas Normalizaci´on de datos 2 horas 3 horas Obtenci´on de datos 55 horas 84 horas Obtenci´on de datos precisos 30 horas 20 horas Procesado de los datos precisos 5 horas 4 horas Generaci´on de dataset 20 horas 60 horas Entrenamiento de modelos 115 horas 52 horas Preprocesado de los datos 15 horas 10 horas Elecci´on de par´ametros de las redes 10 horas 2 horas Entrenamiento de modelos 90 horas 40 horas Evaluaci´on de resultados 18 horas 20 horas An´alisis de las m´etricas 3 horas 2 horas Desnormalizaci´on de resultados 2 horas 1 hora Correcci´on de SGDP4 3 horas 2 horas Evaluaci´on de las correcciones 10 horas 15 horas Documentaci´on 87 horas 82 horas Redacci´on del estudio te´orico 30 horas 35 horas Redacci´on del resto de la memoria 50 horas 40 horas Elaboraci´on de la presentaci´on 7 horas 7 horas Reuniones 14 horas 26 horas Reuni´on de inicio 1.5 horas 1.5 hora Reuniones semanales de media hora 22 sem →11 horas 23 horas Reuni´on de cierre 1.5 hora 1.5 hora Total 450 horas 481 horas Cuadro 2.1: Tareas y duraci´on Las reuniones ser´an con los tutores, Daniel Ayala Hern´andez y Luis Valencia Cabrera, para tareas de evaluaci´on y asesoramiento. Tambi´en se prev´e que en algunas de ellas est´e presente Rafael Ayala Hern´andez, autor del paquete asteRisk, como experto en la materia. 5 Como puede observarse en la tabla, la planificaci´on temporal real se ha desviado de la prevista en un total de 31 horas. Esto ha sido causado, principalmente, por las complicaciones en las tareas de optimizaci´on de los par´ametros orbitales y en las de generaci´on del conjunto de datos. La duraci´on de los sprints scrum que se llevaron a cabo se muestra en la siguiente tabla: Iteraci´on Inicio - Fin Dedicaci´on Art´ıculos y fundamentos matem´aticos 18/11/2022 - 10/12/2022 49 horas Propagaci´on de ´orbitas 11/12/2022 - 09/01/2023 65 horas Redacci´on del estudio te´orico 10/01/2023 - 21/01/2023 35 horas Optimizaci´on de par´ametros 22/01/2023 - 21/02/2023 71 horas Obtenci´on de datos 22/02/2023 - 04/03/2023 31 horas Generaci´on del dataset 05/03/2023 - 02/04/2023 111 horas Entrenamiento y evaluaci´on de modelos 03/04/2023 - 03/05/2023 72 horas Finalizaci´on de la memoria 04/05/2023 - 09/06/2023 40 horas Elaboraci´on de la presentaci´on 10/06/2023 - 16/06/2023 7 horas Cuadro 2.2: Iteraciones Scrum 2.4. Presupuesto Para el c´alculo de costes del proyecto se tendr´an en cuenta varias partidas que ir´an destinadas a cubrir los gastos del personal, la amortizaci´on del equipo inform´atico utilizado y el coste de la electricidad y la conexi´on a Internet. Los costes de personal los vamos a estimar para un desarrollador contratado a jornada completa que perciba una retribuci´on bruta anual de 21.840 euros, lo que supone un coste mensual de 12 euros por hora trabajada. Dentro de los costes del equipo inform´atico hay que considerar la amortizaci´on del equipo empleado durante los 8 meses de duraci´on del proyecto, con un coste inicial de 1245 euros, cuya vida ´util se estima en cuatro a˜nos (48 meses), de los cuales ocho meses ser´an de uso para el proyecto y con un consumo el´ectrico aproximado de 200 W/hora. El coste medio del kilovatio por hora se estima en 0.30 euros y para el suministro de conexi´on a internet se contratar´a una tarifa de 30 euros al mes. Adem´as de los gastos previstos se establece un coste de reservas de riesgos de un 10 por ciento del total para solventar las incidencias que puedan surgir durante la ejecuci´on. De esta manera los costes previstos ascienden a un total de 6221’88e: 6 Partida C´alculo Total Coste de personal 12e/ h x 450 horas 5400 e Amortizaci´on del equipo 1245 e/ 48 meses x 8 meses de proyecto 207’5e Suministro el´ectrico 0’2 kWh x 0’30 ex 450 horas 27 e Conexi´on a internet 30 e/mes / 720 horas/mes x 450 horas 18’75e Total parcial 5656’25 e Reservas (10 %) 565’63 e Total 6221’88 e Cuadro 2.3: Costes previstos Y los costes reales finalmente fueron de: Partida C´alculo Total Coste de personal 12 e/ h x 481 horas 5572 e Amortizaci´on del equipo 1245 e/ 48 meses x 8 meses de proyecto 207’5e Suministro el´ectrico 0’2 kWh x 0’30 ex 481 horas 28’86 e Conexi´on a internet 30 e/mes / 720 horas/mes x 481 horas 20’04 e Gastos finales 5828’4 e Total previsto con reservas 6221’88 e Super´avit en costes 393’48 e Cuadro 2.4: Costes reales y desviaci´on econ´omica Se concluye, por tanto, que se ha producido una desviaci´on de costes de 172’15 euros con respecto a los gastos previstos sin la reserva, lo que representa un 3’04 % de desviaci´on entre dichos costes. Gracias a las reservas, se han podido cubrir dichos gastos y finalmente hay un super´avit de 393’48 euros, lo que representa el 6’32 % del presupuesto total. 7 3. Estudio Te´orico En este cap´ıtulo se desarrolla un estudio te´orico en el que se introduce la teor´ıa de ´orbitas y las bases de los principales modelos de propagaci´on de trayectorias. Dicho estudio est´a basado en [1], [9] y [15]. 3.1. Conceptos previos En esta secci´on se van a recordar y definir conceptos importantes de cara a la comprensi´on del estudio te´orico y desarrollo del proyecto. 3.1.1. Geometr´ıa de la elipse La comprensi´on de esta geometr´ıa resulta fundamental en este trabajo debido a que las ´orbitas que describen los sat´elites alrededor de la Tierra son elipses. Una elipse es el lugar geom´etrico de los puntos de un plano, tales que la suma de sus distancias a dos puntos fijos, llamados focos, siempre es constante. A esta longitud constante se le denomina eje mayor. La mitad de esta longitud se denomina semieje mayor y se denota por a. Al segmento de recta perpendicular al eje mayor cuyos extremos se localizan sobre la elipse se denomina eje menor. Su valor es necesario como dato para poder obtener la ecuaci´on de la elipse. Se denota por ba la mitad de dicho eje, conocido como semieje menor. La excentricidad,e, es una cantidad constante para cada elipse, y da una medida de c´omo de “achatada” es dicha elipse. Se calcula dividiendo la semidistancia focal (de foco a centro), c, entre la longitud del semieje mayor: e=c a=√a2−b2 a Algunos autores prefieren utilizar el achatamiento, en lugar de la excentricidad, para describir la forma. Su f´ormula es: f=a−b a La relaci´on entre ambos par´ametros viene dada por la expresi´on: e2= 2f−f2 Tambi´en se puede definir la elipse como el lugar geom´etrico de los puntos de un plano para los cuales el cociente entre sus distancias a un punto fijo, foco, y a una 8 recta dada, directriz, permanece constante y es igual a la excentricidad de la misma. Cada foco tiene asociada una directriz. La ecuaci´on general de la elipse es: Ax2+By2+Cx +Dy +F= 0 con A=By del mismo signo. Figura 3.1: Esquema de una elipse 3.1.2. Algunas definiciones A continuaci´on se definen ciertos conceptos cuyo conocimiento es necesario para una comprensi´on adecuada de este estudio. ´ Epoca: instante de tiempo que indica un momento en el que se conocen la posici´on y la velocidad del sat´elite. Perturbaci´on: modificaci´on que experimenta el movimiento de un astro a lo largo de su ´orbita como consecuencia de la atracci´on ejercida por los astros pr´oximos. ´ Orbita osculante: ´orbita de Kepler gravitacional (es decir, una el´ıptica u otra c´onica) que tendr´ıa el objeto en el espacio alrededor de su cuerpo central si no hubiera perturbaciones. Anomal´ıa media: ´angulo que forma con el eje de la elipse un sat´elite ficticio que gira con movimiento uniforme sobre una circunferencia cuyo di´ametro coincide con el eje principal de la elipse y llamada circunferencia principal. Anomal´ıa exc´entrica: ´angulo medido desde el centro de la elipse, que forma la proyecci´on del sat´elite sobre la circunferencia principal, y el eje de la elipse. 9 Libraci´on lunar: conjunto de movimientos de oscilaci´on que presenta el disco de la Luna con respecto a un observador ubicado en la Tierra. Ascensi´on recta,α: coordenada astron´omica utilizada para localizar los astros sobre la esfera celeste. Es equivalente a la longitud terrestre (coordenada geogr´afica). Declinaci´on,δ: ´angulo que forma un astro con el ecuador celeste. Es equivalente a la latitud. Acimut: ´angulo de una direcci´on medido en sentido horario a partir del norte geogr´afico. Apoastro, apoapsis, apo´apside, o apocentro: punto de una ´orbita el´ıptica m´as alejado de su centro gravitatorio. En el caso del Sol se conoce como afelio, de la Tierra como apogeo y de la Luna como aposelenio. Periastro, periapsis, peri´apside, o pericentro: punto de una ´orbita el´ıptica m´as cercano a su centro gravitatorio. En el caso del Sol se conoce como perihelio, de la Tierra como perigeo y de la Luna como periselenio. 3.1.3. Leyes de Kepler Estas leyes describen el movimiento de los cuerpos celestes en el espacio. Aunque originalmente formuladas en torno a los planetas y el Sol, sirven para cualquier cuerpo que orbite alrededor de otro. Primera Ley de Kepler (1609) ((Todos los planetas se desplazan alrededor del Sol describiendo ´orbitas el´ıpticas. El Sol se encuentra en uno de los focos de la elipse.)) Segunda Ley de Kepler (1609) ((El radio vector que une un planeta y el Sol recorre ´areas iguales en tiempos iguales.)) Tercera Ley de Kepler (1619) ((Para cualquier planeta, el cuadrado de su per´ıodo orbital es directamente proporcional al cubo de la longitud del semieje mayor de su ´orbita el´ıptica.)) T2 a3=4π2 GM donde Ges la constante de Gravitaci´on Universal y Mes la masa del planeta. 10 3.1.4. Ecuaci´on de Kepler La consecuencia de la segunda ley de Kepler es que el cuerpo no se mueve con un movimiento uniforme, lo cu´al dificulta el c´alculo de la posici´on del planeta dentro de su ´orbita. Por ello aparece la llamada ecuaci´on de Kepler, que nos permite determinar la posici´on exacta en un instante de tiempo determinado, relacionando la anomal´ıa media, M, de c´alculo inmediato, con la anomal´ıa exc´entrica, E, cuyo c´alculo no es tan directo y que es necesario en los modelos de propagaci´on SGP4 y SDP4. A continuaci´on se deduce dicha ecuaci´on. Consideremos la ´orbita definida por una elipse de semieje mayor ay menor b. La masa central, el Sol (F), en torno a la cual orbita la Tierra (Q) se sit´ua, de acuerdo con la primera ley de Kepler, en uno de los focos de la elipse, F. Consideremos tambi´en la circunferencia circunscrita a dicha elipse, cuyo radio es ay tal que, si el planeta recorre la elipse completa en un per´ıodo T, un planeta ficticio recorrer´a la circunferencia completa en el mismo tiempo. Figura 3.2: Situaci´on de una ´orbita Sea el perihelio, P, el punto de partida del movimiento el´ıptico de la Tierra. El ´area barrida en un instante t ser´a la superficie delimitada por PFQ. Vamos a relacionar este ´area con la barrida desde Fpor un planeta ficticio de referencia, Q′, partiendo de Pen su movimiento circular uniforme. Dicha ´area ser´a la limitada por P CS. Por la proporcionalidad de ´areas entre elipse y circunferencia sabemos que: PQ1Q′=PQ1Qa b y FQ1Q′=F Q1Qa b 11 Por tanto, PQ1Q−F Q1Q′= (PQ1Q−F Q1Q)a b y, as´ı, PFQ′=PFQa b Si el per´ıodo orbital es T, en un tiempo t el ´area barrida en la elipse ser´a: PFQ =t Tπab, por lo que en la circunferencia ser´a PFQ′=t Tπa2(I) Por otra parte, PFQ′=PCQ′−FCQ′ y PCQ′=E 2ππa2=E 2a2, donde Eest´a en radianes. Adem´as, FCQ′=1 2ac sen E=1 2a2esen E Luego PFQ′=E 2a2−1 2a2esen E(II) Igualando (I) y (II) llegamos a t T2π=E−esen E y, como M=t T2πse tiene M=E−esen E, que es la denominada ecuaci´on de Kepler. 3.1.5. Constelaciones de sat´elites Hay diversas formas de agrupar los sat´elites, en este trabajo nos interesa mediante constelaciones, que son conjuntos de sat´elites que trabajan como un ´unico sistema, proporcionando una cobertura global o cuasiglobal en todo momento. Entre ellas destacamos: GPS (Global Positioning System), de Estados Unidos. 12 GLONASS, desarrollado por la Uni´on Sovi´etica y administrado hoy en d´ıa por Rusia. EGNOS (European Geostationary Navigation Overlay Service), de la Uni´on Europea. WASS (Wide Area Augmentation System) de Estados Unidos Galileo, de la Uni´on Europea. BeiDou, de China. La importancia de la clasificaci´on seg´un su constelaci´on se comprender´a mejor posteriormente, cuando se expongan los formatos de intercambio de archivos. 3.2. Sistemas de referencia de coordenadas Antes de abordar el estudio de los modelos de propagaci´on conviene definir alguno de los sistemas de coordenadas que van a ser utilizados. Nos centraremos en sistemas de referencia geoc´entricos, puesto que los sat´elites que vamos a estudiar orbitan alrededor de la Tierra, inerciales, es decir, sin aceleraciones del sistema. Trabajaremos con sistemas rectangulares o cartesianos, para los que adem´as del origen hay que indicar el plano fundamental y su vector normal o un sistema ortonormal de vectores. La Tierra y su ´orbita alrededor del Sol forman la base de los sistemas de coordenadas celestes que se van a considerar en el trabajo. Algunos de los m´as destacados son: TEME (True Equator, Mean Equinox): es un sistema geoc´entrico inercial. Resulta complicado encontrar una definici´on exacta de este sistema en la literatura, pero podemos definirlo como aquel en el cual el eje Z coincide con el eje de rotaci´on de la Tierra, el eje X depende del equinocio vernal medio y el eje Y es el producto de los ejes X y Z. Por ejemplo, NORAD lo utiliza en sus archivos TLE. Sistema de Referencia Geoc´entrico Celestial,GCRF: es un marco de coordenadas inercial centrado en la Tierra, donde el origen se coloca en el centro de masa de la Tierra y el marco de coordenadas se fija con respecto a las estrellas, no respecto a la superficie del planeta. El eje X est´a alineado con el equinocio medio de la Tierra a las 12:00 hora terrestre del 1 de enero del a˜no 2000, y el eje Z est´a alineado con el eje de rotaci´on de la Tierra. Es casi equivalente al marco de referencia J2000 (tambi´en llamado EME2000), y en algunos contextos tambi´en se le conoce como marco ICRF (aunque en su definici´on estricta, el origen de coordenadas de ´este ´ultimo se ubica en el baricentro del Sistema Solar). Sistema de Referencia Terrestre Internacional,ITRF: al igual que el GCRF, es un marco de coordenadas inercial centrado en la Tierra y con el origen en el centro de masa de la Tierra. Sin embargo, en este caso el marco gira con respecto a las estrellas, para permanecer fijo con respecto a la superficie de la Tierra a medida que ´esta gira. El eje Z se extiende a lo largo del norte verdadero definido por el polo de referencia IERS, y el eje X se extiende hacia la intersecci´on 13 entre el ecuador y el meridiano de Greenwich en cualquier momento. Tambi´en se conoce como ECEF (Earth-Centered, Earth-Fixed). Cuando se est´en realizando operaciones con posiciones y velocidades hay que poner especial atenci´on en que el sistema de referencia sea el adecuado. 3.3. Sistemas de tiempo El prop´osito del tiempo es definir con precisi´on el momento de un determinado fen´omeno. Por lo tanto, los sistemas de tiempo son el equivalente a los sistemas de referencia pero en el aspecto temporal: fijado un instante determinado, ya podemos medir tiempos relativos a partir de ´el. En el libro ”Principios Matem´aticos de la filosof´ıa natural”(1687), Isaac Newton describe dos tipos de tiempo: absoluto, que no puede ser alterado, y relativo, dependiente de un punto prefijado en la l´ınea temporal. El sistema m´as utilizado es el Tiempo Universal Coordinado, UTC. Surgi´o en 1962, fruto de la cooperaci´on internacional en un esfuerzo por crear un formato de tiempo est´andar y consistente, ya que los que se utilizaban sufr´ıan retrasos debido a las variaciones de la rotaci´on terrestre. El sistema UTC fue finalmente adoptado en 1972. Figura 3.3: Zonas UTC Los sat´elites suelen utilizar el sistema de tiempo GPS que, como su nombre indica, se introdujo con el sistema de navegaci´on que tiene el mismo nombre. La ´epoca (punto de referencia) GPS es el 6 de enero de 1980 a las 0 horas UTC. Debido a la rotaci´on de la Tierra este sistema difiere ligeramente del UTC, por lo que debe ser reajustado diariamente. 14 Ahora pasamos al modelo. Primero recuperamos el movimiento medio (n′′ 0) y el semieje mayor (a′′ 0) originales: a1=ke n02 3 δ1=3 2 k2 a2 1 3 cos2i0−1 (1 −e2 0)3 2 a0=a11−1 3δ1−δ2 1−134 81 δ3 1 δ0=3 2 k2 a2 0 3 cos2i0−1 (1 −e2 0)3 2 n′′ 0=n0 1 + δ0 a′′ 0=a0 1−δ0 Si el perigeo de la ´orbita se encuentra entre 98 y 156 kil´ometros, el valor de la constante s debe reemplazarse por: s∗=a′′ 0(1 −e0)−s+aE Para perigeo inferior a 98 kil´ometros, el valor de s se debe cambiar por: s∗=20 XKMP ER +aE Ahora pasamos a calcular algunas constantes, utilizando los valores apropiados de s: θ= cos i0 ξ=1 a′′ 0−s β0= (1 −e2 0)1 2 η=a′′ 0e0ξ C2= (q0−s)4ξ4n′′ 0(1 −η2)−7 2a′′ 01 + 3 2η2+ 4e0η+e0η3 +3 2 k2ξ (1 −η2)−1 2+3 2θ2(8 + 24η2+ 3η4) C1=B∗C2 C4= 2n′′ 0(q0−s)4ξ4a′′ 0β2(1 −η2)−7 22η(1 + e0η)1 2e0+1 2η3−2k2ξ a′′ 0(1 −η2)× 3(1 −3θ2)1 + 3 2η2−2e0η3+3 4(1 −θ2)(2η2−e0η−e0η3) cos 2ω0 ˙ M=1 + 3k2(−1+3θ2) 2a′′ 2 0β3 0 +3k2 2(13 −78θ2+ 137θ4) 16a′′ 4 0β7 0n′′ 0 21 ˙ω=−3k2(−1−5θ2) 2a′′ 2 0β4 0 +3k2 2(7 −114θ2+ 395θ4) 16a′′ 4 0β8 0 +5k2(3 −36θ2+ 49θ4) 4a′′ 4 0β8 0n′′ 0 ˙ Ω1=−3k2θ a′′ 4 0β4 0 n′′ 0 ˙ Ω = ˙ Ω1+3k2 2(4θ−19θ3 2a′′ 4 0β8 0 +5k4θ(3 −7θ2) 2a′′ 4 0β8 0n′′ 0 En este momento se realiza la inicializaci´on de los par´ametros que sirven para calcular las ecuaciones de perturbaci´on del espacio profundo. Los efectos de la gravedad se incluyen en: MDF =M0+˙ M(t−t0) ωDF =ω0+ ˙ω(t−t0) ΩDF = Ω0+˙ Ω(t−t0) El efecto del rozamiento en la longitud del nodo ascendente viene determinado por: Ω = ΩDF −21 2 n′′ 0k2θ a′′ 2 0β2 0 C1(t−t0)2 A continuaci´on habr´ıa que a˜nadir los efectos del espacio profundo y de la resonancia de largo periodo a los seis elementos orbitales cl´asicos: a=aDS[1 −C1(t−t0)]2 e=eDS −B∗C4(t−t0) L=MDS +ωDS + ΩDS +n′′ 03 2C1(t−t0)2 donde aDS,eDS,MDS,ωDS y ΩDS son los valores de n0,e0,MDF ,ωDF y ΩDF despu´es de aplicar las perturbaciones de a resonancia y el espacio profundo. Ahora hay que incluir las perturbaciones de espacio profundo debidas al Sol y la Luna. A partir de aqu´ı, denotaremos por n,e,I,ω, Ω y Mal movimiento medio, la excentricidad, la inclinaci´on, el argumento de perigeo, la longitud del nodo ascendente y la anomal´ıa media, respectivamente, despu´es de a˜nadir los periodos solares y lunares. axN =ecos ω β=√1−e2 LL=A3,0sin i0 8k2aβ2ecos ω3+5θ 1 + θ ayNL =A3,0sin i0 4k2aβ2 LT=L+LL ayN =esin ω+ayNL 22 Resolvemos la ecuaci´on de Kepler para E+ωdefiniendo U=LT−Ω y aplicando el m´etodo de aproximaciones sucesivas: (E+ω)i+1 = (E+ω)i+ ∆(E+ω)i donde ∆(E+ω)i+1 =U−ayN cos (E+ω)i+axN sin (E+ω)i−(E+ω)i −ayN sin (E+ω)i−axN cos (E+ω)i+ 1 y (E+ω)1=U Las siguientes ecuaciones se utilizan para calcular valores que necesitaremos con posterioridad. ecos E=axN cos (E+ω+ayN sin (E+ω) esin E=axN sin (E+ω+ayN cos (E+ω) eL= (a2 xN +a2 yN )1 2 pL=a(1 −e2 L) r=a(1 −ecos E) ˙r=ke √a resin E r˙ f=ke √pL r cos u=a r"cos (E+ω)−axN +ayN (esin E) 1 + p1−e2 L# sin u=a r"sin (E+ω)−ayN +axN (esin E) 1 + p1−e2 L# u= tan−1sin u cos u ∆r=k2 2pL (1 −θ2) cos 2u ∆u=−k2 4pL2(7θ2−1) sin 2u ∆Ω = 3k2θ 2pL2sin 2u ∆i=3k2θ 2pL2sin i0cos 2u 23 ∆ ˙r=−k2n pL (1 −θ2) sin 2u ∆r˙ f=k2n pL(1 −θ2) cos 2u−3 2(1 −3θ2) A˜nadimos estos incrementos para obtener las cantidades osculantes: rk=r"1−3 2k2p1−e2 L pL2(3θ2−1)#+ ∆r uk=u+ ∆u Ωk= Ω + ∆Ω ik=I+ ∆i ˙rk= ˙r+ ∆ ˙r r˙ fk=r˙ f+ ∆r˙ f Los vectores unitarios de orientaci´on se calculan con  U= Msin uk+ Ncos uk  V= Mcos uk+ Nsin uk donde  M= −sin Ωkcos ik cos Ωkcos ik sin ik  y  N=  cos Ωk sin Ωk 0  Finalmente, la posici´on r y la velocidad ˙ r se obtienen como r =rk U y˙ r = ˙rk U+r˙ fk V 3.7. Modelos de propagaci´on de alta precisi´on Johann Franz Encke (1791-1865) present´o un m´etodo de c´alculo de ´orbitas para el caso en el que las perturbaciones fueran peque˜nas con respecto al movimiento de los cuerpos. Dicho m´etodo se basaba en integrar num´ericamente las diferencias entre 24 la ´orbita osculante y la ´orbita real. De este modo se logr´o suficiente precisi´on para las limitadas capacidades computacionales de la ´epoca. Philip Herbert Cowell (1870–1949) utiliz´o una t´ecnica num´erica para determinar la ´orbita del octavo sat´elite de J´upiter. Con el aumento de la potencia computacional, su t´ecnica gan´o cada vez m´as popularidad hasta el punto de que hoy en d´ıa, en astrodin´amica, se utiliza la formulaci´on de Cowell para establecer las ecuaciones del movimiento, para posteriormente realizar integraci´on num´erica. 3.7.1. Formulaci´on de Encke Gracias al aumento de la capacidad computacional, este m´etodo est´a en desuso, pero debido a su importancia hist´orica merece una menci´on. No se basa en integrar todas las fuerzas que act´uan sobre el sat´elite, sino en integrar la diferencia entre la aceleraci´on debida al astro que est´a en el foco y la aceleraci´on perturbada, siempre y cuando esta sea peque˜na. As´ı se consigue una precisi´on buena sin mucho coste. Figura 3.6: Formulaci´on de Encke La aceleraci´on y la aceleraci´on perturbada son: ¨ r =−µ r3r ¨ rp=−µ r3 p rp+ap donde µ=GM, con G≈6,67 ×10−11 N·m2 kg es la constante de gravitaci´on universal Mes la masa de la Tierra 25 Por tanto, la diferencia entre las aceleraciones es: δ¨ r =ap+µ r31−r3 r3 prp−δr 3.7.2. Formulaci´on de Cowell La idea detr´as de este m´etodo es simple pero muy potente, gracias al aumento de la capacidad computacional. Se basa en incorporar aceleraciones a la ecuaci´on del movimiento de dos cuerpos para obtener una ecuaci´on m´as precisa sobre el movimiento: ¨ r =−µ r3r +ap donde apes la aceleraci´on total causada por todas las fuerzas que act´uan sobre el sat´elite. La forma concreta de apdepender´a del n´umero de perturbaciones que estemos considerando. Ahora el problema que se nos presenta es la obtenci´on de los datos necesarios para implementar el modelo que queramos. Por ejemplo, no podemos determinar con exactitud los coeficientes gravitacionales de Plut´on (debido a que no conocemos con exactitud su radio ecuatorial), luego si los incluimos en el modelo estaremos realizando peque˜nas aproximaciones. La ecuaci´on anterior se conoce como formulaci´on de Cowell. Conviene distinguir entre la formulaci´on de Cowell y el m´etodo de Cowell. Este ´ultimo consiste en aplicar el m´etodo de las diferencias finitas para integrar la ecuaci´on. Normalmente se reformulan las tres ecuaciones de segundo orden del movimiento en un sistema de seis ecuaciones de primer orden, de modo que el problema pueda ser abordado con una gama m´as amplia de m´etodos de integraci´on num´erica. Dicho sistema se conoce como una variaci´on de la formulaci´on de Cowell:  X=r v  ˙  X="v −µr r3+ap# 3.7.3. HPOP El m´etodo de propagaci´on de alta precisi´on HPOP se basa en la formulaci´on de Cowell, teniendo en cuenta una multitud de fuerzas que act´uan sobre el sat´elite. Dicha ecuaci´on quedar´ıa del siguiente modo: ¨ r =−µ r3r +X cc∈CC acc +adrag +asrp 26 donde CC son los cuerpos celestes del Sistema Solar distintos a la Tierra, es decir, el Sol, la Luna y todos los planetas menos la Tierra. A continuaci´on se explican estas perturbaciones adicionales, no sin antes mencionar que las coordenadas y velocidades ingresadas y calculadas con este m´etodo est´an en el marco de referencia GCRF. Campos gravitatorios Se consideran los efectos de los campos gravitatorios del Sol, la Luna y el resto de planetas del Sistema Solar. Esta aceleraci´on la podemos expresar como: acc =−GM r2 r ∥r∥ donde Ges la constante de Gravitaci´on Universal Mes la masa del cuerpo celeste r es el vector posici´on del sat´elite respecto del cuerpo celeste Rozamiento atmosf´erico La ecuaci´on b´asica del rozamiento aerodin´amico es: adrag =−1 2 cDA mρv2 rel vrel ∥vrel∥ donde cDes el coeficiente de rozamiento ρes la densidad atmosf´erica Aes el ´area de la secci´on transversal, es decir, el ´area de la secci´on normal al vector velocidad del sat´elite mes la masa del sat´elite vrel es la velocidad relativa a la atm´osfera Cabe mencionar que normalmente m cDAse conoce como coeficiente bal´ıstico, BC. Si queremos calcular la velocidad relativa a la atm´osfera, sin tener en cuenta el movimiento del viento, podemos utilizar la siguiente expresi´on: 27 vrel =dr dt −ω⊕×r =           dx dt +ω⊕y dy dt −ω⊕x dz dt           Si tenemos en cuenta el movimiento del viento la expresi´on quedar´ıa: vrel =           dx dt +ω⊕y+vw[−cos(α)sen(δ)cos(βw)−sen(α)sen(βw)] dy dt −ω⊕x+vw[−sen(α)sen(δ)cos(βw)−cos(α)sen(βw)] dz dt + +vw[cos(δ)cos(βw)]           donde αes la ascensi´on recta δes la declinaci´on vwes la velocidad del viento βwes el acimut Presi´on de la radiaci´on solar La aceleraci´on debida a la radiaci´on solar se puede expresar como: asrp =−psrpcRA⊙ m rsat⊙ ∥rsat⊙∥ donde psrp es la presi´on solar cRes la reflectividad A⊙es el ´area de secci´on transversal expuesta al Sol rsat⊙es el vector velocidad respecto al Sol 28 Mareas Podemos distinguir entre mareas terrestres, que pueden ser deformaciones de la Tierra causadas principalmente por la atracci´on gravitatoria de la Luna y el Sol o por el efecto centr´ıfugo de la rotaci´on, y mareas oce´anicas, especialmente influenciadas por la Luna. Las expresiones matem´aticas de estas perturbaciones se pueden expresar mediante polinomios de Legendre y son muy complejas, por lo que se omiten en el presente trabajo. 3.7.4. Resoluci´on num´erica de ecuaciones diferenciales Los m´etodos de resoluci´on num´erica de ecuaciones diferenciales forman un conjun- to de herramientas matem´aticas que nos permiten encontrar aproximaciones num´ericas a las soluciones de ecuaciones diferenciales. Tambi´en se conoce como integraci´on num´erica aunque este nombre puede confundirse con el c´alculo aproximado de integrales definidas. Su importancia radica en que la mayor´ıa de las veces no se puede resolver anal´ıticamente las ecuaciones diferenciales. Para profundizar en el tema se puede consultar [4] y [8]. A continuaci´on se presentan el m´etodo de Runge-Kutta de cuarto orden, que es el que utiliza la implementaci´on de HPOP de asteRisk, y el m´etodo de Runge-Kutta- Felhberg, que supone una mejora del anterior. M´etodo de Runge-Kutta (de 4ºorden) Posiblemente las f´ormulas de Runge-Kuta son los m´etodos de integraci´on num´erica m´as conocidos, desarrollados a principios del siglo XX por Carl Runge y Wilhen Kutta, siendo aplicaciones de las series de Taylor. A continuaci´on se describir´a el m´etodo de cuarto orden. Consideremos la ecuaci´on diferencial de primer orden y′=f(x, y) con la condici´on inicial y(x0) = y0 y siendo Consideremos el siguiente problema de valores iniciales: (y′=f(x, y) y(x0) = y0 donde f∈C5, es decir, la quinta derivada de fexiste y es continua y elijamos una anchura de paso hy sean xi=x0+ih eyi=y(xi), para i= 0,1,2.... Consideremos: 29                          k(i) 1=f(xi, yi) k(i) 2=f xi+h 2, yi+k(i) 1 2! k(i) 3=f xi+h 2, yi+k(i) 2 2! k(i) 4=fxi+h, yi+k(i) 3 As´ı, el procedimiento de Runge-Kutta de cuarto orden es yi+1 =yi+1 6k(i) 1+ 2k(i) 2+ 2k(i) 3+k(i) 4 Figura 3.7: Representaci´on gr´afica del m´etodo de Runge-Kutta de cuarto orden M´etodo de Runge-Kutta-Felhberg Tambi´en conocido como m´etodo de Runge-Kutta encajado om´etodo de Felhberg, emplea un tama˜no de paso variable. Supone una buena mejora del algoritmo ajustando el valor del paso para mantener peque˜nos los errores locales de truncamiento. Consiste en aplicar los m´etodos de cuarto y quinto orden de Runge-Kutta y compararlos. Este m´etodo es muy bueno para utilizar en la cerca del apoastro, pero no se recomienda su uso cerca del periastro, donde el sat´elite se mueve m´as r´apido. Se consideran los siguientes valores: 30 El repositorio deb´ıa contener archivos RINEX versi´on 2.10, que es la soportada por el paquete asteRisk, por lo que el escogido fue el del Instituto de Estad´ıstica y Cartograf´ıa de la Junta de Andaluc´ıa [3]. Concretamente, se utilizaron los datos procedentes del receptor de la estaci´on de C´ordoba, situada en la Consejer´ıa de Fomento y Vivienda. Dicho repositorio dispone de datos sobre sat´elites GPS y GLONASS. Para descargar los archivos se utiliz´o web scrapping. Posteriormente se realiz´o su proceso y an´alisis, para determinar cu´al era el sat´elite con m´as mensajes RINEX disponibles y generar el dataset con sus datos. El procesado se ha realizado con ayuda de la funci´on de lectura RINEX incluida en el paquete asteRisk. Para la elecci´on del sat´elite se procedi´o a la lectura de todos los archivos, contando el n´umero de mensajes RINEX que hab´ıa de cada sat´elite, de modo que al terminar se pod´ıa obtener sin dificultad cu´al era el sat´elite del que m´as informaci´on hab´ıa. Cabe mencionar que el tratamiento de los GPS y los GLONASS hab´ıa que hacerlo por separado, ya que los campos de la lista devuelta por las respectivas funciones de lectura eran distintos. El sat´elite con m´as informaci´on disponible result´o ser el Kosmos 2514, de la constelaci´on GLONASS y con c´odigo identificativo 17. Para almacenar sus datos en un dataframe se realiz´o de nuevo la lectura de todos los archivos y se fue guardando la posici´on y velocidad, separadas en coordenadas, y la ´epoca de los mensajes correspondientes a este sat´elite, convirtiendo previamente la posici´on y velocidad al marco de referencia GCRF (puesto que estaban en el ITRF). As´ı, dicho conjunto de datos ten´ıa la siguiente estructura: ´ Epoca Posici´on Velocidad Cuadro 5.1: Datos obtenidos de los ficheros RINEX del sat´elite Una vez tuvimos el dataframe con las posiciones, velocidades y ´epocas pasamos a la generaci´on del dataset para entrenar los modelos. La generaci´on se realiz´o en dos pasos. En primer lugar, se guardaron los datos de la fecha inicial, la posici´on y velocidad iniciales reales, la fecha de la predicci´on, la posici´on y velocidad reales en esta segunda fecha y la posici´on y velocidad predichas en ese instante con los par´ametros optimizados. En el segundo paso, se a˜nadieron las predicciones con los par´ametros sin optimizar. La elecci´on de las fechas base se realiz´o tomando de manera aleatoria 5000 fechas entre las del dataframe del sat´elite. En la optimizaci´on de los elementos orbitales y la elecci´on de fechas para las predicciones se mantuvo la idea inicial, con dos peque˜nas modificaciones, motivadas por la posible falta de datos en una fecha concreta, ya que no sabemos si disponemos de ella exactamente. Estos dos cambios fueron: La optimizaci´on de los par´ametros orbitales en cada punto se realiz´o en base a los 5 anteriores, ´el mismo y los 5 posteriores, ordenados por fecha. En el tiempo a predecir se permiti´o una desviaci´on del 30 %, es decir, para predecir a 100 minutos, se utiliz´o la fecha del conjunto m´as cercana, siempre y cuando ´esta distara un m´aximo de 130 minutos de la correspondiente. De este modo, el dataset tom´o esta estructura, donde si aparece “optim” quiere decir que se han empleado los par´ametros optimizados: 37 ´ Epoca inicial Posici´on inicial Velocidad inicial ´ Epoca de la predicci´on Posici´on predicha optim Velocidad predicha optim Cuadro 5.2: Dataset con datos reales y predicciones SGDP4 con par´ametros optimizados Una vez tuvimos este dataset faltaba a˜nadir las predicciones de SGDP4 con los par´ametros sin optimizar. Para ello se iba recorriendo las filas de dicho conjunto de datos y con los valores de la ´epoca, posici´on y velocidad iniciales, y la ´epoca de la predicci´on realizamos la predicci´on en cuesti´on y la a˜nadimos al dataset, adquiriendo la siguiente estructura: ´ Epoca inicial Posici´on inicial Velocidad inicial ´ Epoca de la predicci´on Posici´on predicha optim Velocidad predicha optim - Posici´on predicha Velocidad predicha Cuadro 5.3: Dataset con posiciones reales y predicciones SGPD4 con par´ametros optimizados y sin optimizar 5.3. Creaci´on de los modelos Para la creaci´on y entrenamiento de modelos se opta por TensorFlow y Keras con R, por lo que hay que utilizar la versi´on de estas librer´ıas en dicho ecosistema. Adem´as, se considera el uso de CUDA y cuDNN para acelerar el proceso de entrenamiento. Para la separaci´on de los datos en conjuntos de entrenamiento y test se utiliz´o la librer´ıa rsample. Antes de comenzar con el entrenamiento de los modelos debemos normalizar los datos, es decir, reducir su rango y que as´ı est´en aproximadamente en la misma escala, lo que ayudar´a a la red a realizar mejores predicciones. La normalizaci´on elegida fue el escalado robusto (Robust Scaler), que consiste en sustraer la mediana y escalar los datos de acuerdo al rango intercuart´ılico (IQR). Si lo queremos expresar en una f´ormula quedar´ıa: xi−Q2(x) Q3(x)−Q1(x) La ventaja del escalado robusto radica en que los valores que utiliza para escalar (mediana y rango intercuart´ılico) no se ven afectados por los valores at´ıpicos (extremos), luego mantiene la distancia relativa entre los valores pr´acticamente intacta. Notar que los par´ametros de este m´etodo de normalizaci´on, es decir, la mediana y el rango intercuart´ılico, se calculan sobre el conjunto de entrenamiento y son los que se le aplica al conjunto de test y a cualquier otro dato de entrada con el que se vayan a realizar predicciones. La estructura de la red neuronal implementada y entrenada se compone de una capa de entrada, con tantas neuronas como datos de entrada haya, 3 capas densas 38 con activaci´on relu y 256, 128 y 64 neuronas, respectivamente, y la capa de salida con 3 neuronas, una por cada coordenada (X, Y, Z), sin funci´on de activaci´on (i.e. la identidad), ya que se quiere realizar regresi´on para estimar el error de la predicci´on para cada coordenada. Figura 5.1: Estructura de las redes neuronales implementadas Como funci´on de p´erdida se eligi´o el error cuadr´atico medio, MSE. La t´ecnica de optimizaci´on elegida fue ADAM (Momento Adaptativo) y la m´etrica fue el error absoluto medio, MAE, para tener una idea de cu´anto nos desviamos del valor real en su propia escala. Tras varias pruebas con algunas combinaciones, el tama˜no de batch elegido fue de 128 y el n´umero de epochs del entrenamiento de 100. La divisi´on en conjuntos de entrenamiento y prueba se realiz´o con estratificaci´on, considerando el tiempo a predecir, es decir, la diferencia entre la fecha de la predicci´on y la fecha base. 5.4. Evaluaci´on de resultados Con la m´etrica empleada MAE podemos tener cierta idea del rendimiento del modelo, pero no tiene una interpretaci´on f´acil debido a que los datos se encuentran normalizados. Por tanto, se decide desnormalizar la salida de la red, para obtener una predicci´on del error en la escala original y, posteriormente, sumar ese error a la predicci´on de SGDP4 para obtener su correcci´on. Tras obtener la predicci´on de la posici´on corregida se procede a calcular el error medio con respecto a las posiciones reales y compararlo con el de las predicciones sin corregir, para cuantificar realmente si mejora la predicci´on de la posici´on y cu´anto lo hace. 39 5.5. Ampliaci´on de la experimentaci´on Despu´es de obtener los primeros modelos se procedi´o a˜nadir datos sobre los astros del sistema solar al dataset, ya que sus campos gravitatorios y movimientos oscilatorios (como la libraci´on lunar) tienen efecto sobre los sat´elites. La obtenci´on de estos datos se realiz´o utilizando el modelo JPL DE440, proporcionado por el Jet Propulsion Laboratory de la NASA e incluido como una funci´on en el paquete asteRisk bajo el nombre de JPLephemerides. Previamente a la adici´on al conjunto de datos de los sat´elites, se cre´o otro dataset auxiliar para los astros, con la posici´on y velocidad de todos los astros del Sistema Solar salvo la Tierra (puesto que al estar en el marco GCRF valdr´ıan 0) y la libraci´on lunar y sus derivadas primera y segunda, en todas las fechas que se encontraban entre las fechas base del dataset de los sat´elites. Las posiciones y velocidades separadas en sus coordenadas de los ejes X, Y y Z, y la libraci´on lunar en sus componentes ϕ,θy ψ. De este modo, la estructura de este dataset qued´o de la siguiente forma: ´ Epoca Posici´on del Sol Velocidad del Sol - Posici´on de Mercurio Velocidad de Mercurio - Posici´on de Venus Velocidad de Venus - Posici´on de Marte Velocidad de Marte - Posici´on de J´upiter Velocidad de J´upiter - Posici´on de Saturno Velocidad de Saturno - Posici´on de Urano Velocidad de Urano - Posici´on de Neptuno Velocidad de Neptuno - Posici´on de Plut´on Velocidad de Plut´on Libraci´on Lunar d(Libraci´on Lunar) dd(Libraci´on Lunar) Cuadro 5.4: Dataset con datos de los astros Una vez generado este dataset se procede a unirlo con el de las posiciones de sat´elites, para as´ı obtener un conjunto de datos compuesto por: una ´epoca inicial, posici´on y velocidad del sat´elite; los datos mencionados anteriormente de los astros en dicha ´epoca, la ´epoca en la que realizar la predicci´on y la posici´on y velocidad predichas por SGDP4, tanto con los par´ametros optimizados como con los par´ametros sin optimizar en dicha ´epoca. Las fechas mencionadas estaban guardadas en formato cadena de texto (para que sea legible en caso de que haya que depurar) y en formato UNIX, y el resto de datos guardados por componentes y en el marco de referencia GCRF. Este nuevo dataset nos sirvi´o para generar nuevos modelos que recibieran como entrada lo mismo que los anteriores y, adem´as, los datos de los astros en las fechas iniciales. No se guardaron en el dataset los errores de las predicciones, porque tras realizar pruebas de rendimiento, era m´as r´apido leer el fichero en el que estaba guardado y a˜nadirlos seg´un fuera necesario que almacenarlos en el propio dataset. 40 6. Problemas encontrados En esta secci´on se expondr´an las complicaciones que tuvieron lugar durante el desarrollo del proyecto, as´ı como la soluci´on aportada para cada una de ellas. 6.1. Ajuste de par´ametros El primer problema aparece con la optimizaci´on de los par´ametros orbitales. Espec´ıficamente, se detect´o que el paquete no trataba algunos casos minoritarios, la mayor´ıa de los cuales eran con valores cercanos a la frontera de su rango, en los que se produc´ıa una divisi´on entre 0 que causaba error. Dicho error fue reportado para su futuro arreglo por parte del autor del paquete. Para poder realizar el ajuste adecuadamente se propuso utilizar un manejador de excepciones que devolviera un valor muy alto en la funci´on de error a minimizar en caso de fallo, de modo que nunca se eligiera la combinaci´on de par´ametros en cuesti´on. El segundo problema surge en algunas predicciones que a corto plazo se realizaban sin problema, pero cuando se aumentaba unas horas el tiempo de la predicci´on daba error. Se notific´o tambi´en al autor del paquete para su arreglo. Este fallo apareci´o cuando se recorr´ıa una lista de fechas para realizar la predicci´on en ellas y a˜nadirlas al dataset, por lo que la soluci´on implementada fue utilizar un manejador de excepciones que en caso de fallo, diera como predicci´on NA. 6.2. Representaci´on de ´orbitas Al representar las trayectorias obtenidas por SGDP4 con los par´ametros optimizados y sin optimizar, as´ı como las generadas por HPOP, se observ´o que la trayectoria de los par´ametros optimizados difer´ıa excesivamente de las otras dos, lo cual es signo de alg´un problema. Se baraja que puedan ocurrir tres cosas: La funci´on de optimizar no sea correcta No se est´e haciendo bien el cambio de sistemas de coordenadas El m´etodo de estimaci´on de par´ametros iniciales a partir de una posici´on y velocidad iniciales, que hab´ıa sido a˜nadido recientemente al paquete, no fuera correcto Para poder descubrir cu´al era el problema hubo que realizar una profunda investigaci´on, que incluy´o depurar el m´etodo de estimaci´on de par´ametros y multitud de pruebas con la funci´on de optimizar y con cambios de sistemas de coordenadas para comprobar que las transformaciones entre marcos de referencia eran correctas. Durante este proceso se identific´o un error en un cambio de coordenadas que no se estaba haciendo correctamente; sin embargo, el error global no se deb´ıa a ese fallo. Finalmente se descubri´o que el error proven´ıa del formato de fecha, que por defecto es CET o CEST 41 y no UTC (difieren en 1 ´o 2 horas), por lo que al hacer el cambio de sistema no se estaba pasando adecuadamente la ´epoca. 6.3. Repositorios La b´usqueda de un repositorio con archivos RINEX versi´on 2.10, compatible con las funciones de lectura de asteRisk, no fue f´acil, debido a que la mayor´ıa de los repositorios de las agencias estatales est´an en la versi´on 3 y las que se encontraban de versi´on 2 eran de la 2.11 en adelante. En los repositorios de la NASA s´ı que hab´ıa de la versi´on 2.10, pero su estructura de directorios era compleja y no resultaba sencillo distinguir dicha versi´on de la 2.11. Finalmente se encontr´o el repositorio del Instituto de Estad´ıstica y Cartograf´ıa de la Junta de Andaluc´ıa, donde estaban separados completamente los ficheros de versi´on 2 de los de la 3 y, adem´as, los de la versi´on 2 eran concretamente la 2.10. 6.3.1. Estaci´on de la Universidad de Sevilla Los archivos RINEX que se descargaron la primera vez correspond´ıan a la estaci´on de la Universidad de Sevilla, pero durante el generaci´on del conjunto de datos aparecieron errores inesperados que, tras depurar el c´odigo e investigar, resultaron estar causados por la presencia de observaciones de datos repetidas, es decir, varios mensajes RINEX iguales. Se procedi´o de nuevo a la lectura de los datos del sat´elite, puesto que no deber´ıan aparecer mensajes duplicados en dichos archivos, y result´o que todas las fechas, ´epocas de las efem´erides, hab´ıan cambiado. Tras esto se empez´o un an´alisis m´as profundo y durante la lectura manual de un archivo RINEX apareci´o el problema: ciertos campos de la cabecera, tales como las correcciones temporales, se encontraban vac´ıos, lo que ocasionaba que, durante el lectura que realiza la funci´on del paquete, las variables correspondientes tomaran el valor NULL y el c´alculo de la fecha, al aplicar la correcci´on del bias, no se realizaba correctamente. Despu´es de analizar la situaci´on se decidi´o cambiar la estaci´on de la que se descargaron los datos, escogi´endose finalmente la de C´ordoba, cuyos archivos RINEX s´ı dispon´ıan de las cabeceras completas y no mostraban ning´un problema. 42 7. Implementaci´on En este cap´ıtulo se presentar´an las principales herramientas empleadas en el desarrollo del proyecto, as´ı como las funciones implementadas y los extractos de c´odigo elaborados. Entre dichas funciones se encuentran las que han sido utilizadas en la creaci´on del dataset utilizado para entrenar la red neuronal. Menci´on especial merece la creaci´on de dicho conjunto de datos, puesto que supone una contribuci´on de gran peso, que involucra ajustes de par´ametros complejos y se ha estructurado de un modo que simplifica su tratamiento, encontr´andose todos los datos de magnitudes equiparables en el mismo marco de referencia y con las mismas unidades. 7.1. Herramientas 7.1.1. R R es un lenguaje de programaci´on sumamente potente para realizar distintos c´alculos cient´ıficos, num´ericos y estad´ısticos, as´ı como para crear gr´aficas y figuras de gran calidad. Surge de la mano de Robert Gentleman y Ross Ihaka del Departamento de Estad´ıstica de la Universidad de Aucklanden en 1993 como una reimplementaci´on de software libre del lenguaje S, adicionado con soporte para ´ambito est´atico. Se trata de uno de los lenguajes de programaci´on m´as utilizados en investigaci´on cient´ıfica, siendo adem´as muy popular en los campos de aprendizaje autom´atico, miner´ıa de datos, bioinform´atica o econometr´ıa, entre otros. Todo esto gracias a la posibilidad de utilizar multitud de paquetes con amplias funcionalidades de c´alculo y representaci´on gr´afica. Se distribuye bajo la licencia GNU GPL y est´a disponible para los sistemas operativos Windows, Macintosh, Unix y GNU/Linux. Los paquetes de R m´as utilizados para este proyecto se describen seguidamente. 7.1.2. asteRisk y asteRiskData El paquete de R asteRisk proporciona funcionalidades b´asicas para el c´alculo de posiciones de sat´elites a partir de un vector de estado. Incluye implementaciones de: Los modelos de perturbaciones simplificadas SGP4 y SDP4, incluida una funci´on sgdp4 que elige el m´as adecuado de los dos. Un modelo de alta precisi´on (HPOP). Utilidades para leer archivos TLE y RINEX. Cambio de coordenadas entre sistemas de referencias. 43 Modelos de posiciones de los astros, como el JPL DE440. Muchas de las funcionalidades del paquete requieren de multitud de coeficientes y otros datos, como informaci´on sobre las mareas o tormentas solares, por lo que tambi´en es necesario el paquete asteRiskData. 7.1.3. Dplyr Forma parte del ecosistema de paquetes tidyverse, que trabajan conjuntamente para propiciar una ciencia de datos mejor estructurada y que incluye dplyr, tidyr, readr o ggplot2, entre otros. El paquete dplyr es un versi´on optimizada del paquete plyr. Contiene funciones muy potentes que ser´ıa muy complicado replicar, otra contribuci´on importante es que proporciona una ”gram´atica”para manipular y realizar operaciones con dataframes. Esto es muy ´util, ya que proporciona una abstracci´on que anteriormente no exist´ıa y permite que el c´odigo sea m´as legible. Tambi´en cabe destacar que las funciones de dicho paquete son muy r´apidas, puesto que est´an implementadas con el lenguaje C++. 7.1.4. Plotly Con Plotly se pueden construir gr´aficas de alta calidad e incorporar interactividad en los gr´aficos, mostrar la informaci´on que contiene cada dato, mostrar datos filtrando por alg´un factor, etc. Adem´as, los gr´aficos son ajustables (responsive en ingl´es), es decir, se adaptan a las dimensiones de la ventana en que aparecen. 7.1.5. Reticulate Reticulate proporciona una interfaz para utilizar Python desde R. Este paquete ser´a utilizado internamente por las librer´ıas tensorflow y keras de R. 7.1.6. Python Python es un lenguaje de programaci´on enfocado en la legibilidad del c´odigo. Tiene su origen a finales de los a˜nos ochenta en el Stichting Mathematisch Centrum (CWI) de Pa´ıses Bajos. Su creador fue Guido van Rossum, que se inspir´o en el lenguaje ABC y le puso el nombre debido a su afici´on por los humoristas Monty Python. Dispone de un cat´alogo muy amplio de librer´ıas en multitud de campos, desde la inteligencia artificial hasta la optimizaci´on. A continuaci´on se describen las librer´ıas de Python que van a ser utilizadas en este trabajo. 44 7.1.7. Beautiful Soup Es una librer´ıa concebida para analizar textos HTML, incluso aunque el marcado sea incorrecto. Por lo tanto, es ampliamente utilizada para realizar web scrapping. 7.1.8. TensorFlow Es una liber´ıa para aprendizaje autom´atico. Desarrollada por Google y de c´odgio abierto permite la creaci´on y el entrenamiento de redes neuronales de un modo muy c´omodo. Pese a que no est´a adaptado como tal el paquete a R, s´ı existe un paquete que realiza vincula R con Python y permite utilizar TensorFlow con dicho lenguaje. 7.1.9. CUDA y cuDNN Son dos tecnolog´ıas desarrolladas por NVIDIA para acelerar el procesamiento en paralelo en las tarjetas gr´aficas (GPUs) y mejorar as´ı el rendimiento en aplicaciones de aprendizaje profundo (deep learning.) 7.2. Desarrollo 7.2.1. Instalaci´on de asteRisk y asteRiskData Lo primero que debemos hacer es instalar los paquetes de R asteRisk yasteRisk- Data. Para ello ejecutamos los siguientes comandos: 1install.packages("asteRisk") 2install.packages("asteRiskData", repos="https://rafael-ayala.github.io/ drat/") Extracto de c´odigo 7.1: Instalaci´on asteRisk y asteRiskData Cada vez que se carga asteRiskData aparece por consola un aviso de que se ejecute el comando getLatestSpaceData() para obtener los datos m´as recientes de tormentas solares, mareas, etc. La instrucci´on ser´ıa la siguiente: 1getLatestSpaceData(target="all") Extracto de c´odigo 7.2: Obtenci´on de datos recientes de los astros 7.2.2. C´alculo de elementos cl´asicos orbitales El marco de referencia que se decidi´o utilizar es el GCRF y como el paquete asteRisk no dispone de una funci´on que, a partir de la posici´on y velocidad de un 45 sat´elite en dicho sistema, obtenga los elementos orbitales, debemos implementarla. En dicha funci´on ser´a necesario realizar la conversi´on de GCRF a TEME, que involucra el uso de cuaterniones para aplicar las rotaciones pertinentes, y posteriormente utilizar la funci´on ECItoKOE con la conversi´on de unidades de medida adecuadas. 1get_params_orbitales <- function(dateTime, pos, vel){ 2# Comienza la transformacion a TEME 3# Hace falta la fecha en MJD_UTC. 4date_mjd_utc <- dateTimeToMJD(dateTime, timeSystem = "UTC") 5# Cuaternion con "rotacion necesaria" 6quaternion <- asteRisk:::rotationGCRFtoTEME(date_mjd_utc) 7# Convierte el cuaternion a "Directed Cosine Matrix" 8dcm <- asteRisk:::quaternionToDCM(quaternion) 9# Obtencion del TEME, %* % es la multiplicacion de matrices 10 pos <- as.numeric(dcm %* % pos) 11 vel <- as.numeric(dcm %* % vel) 12 koe <- ECItoKOE(pos, vel) 13 # semiMajorAxisToMeanMotion devuelve la rotacion en vueltas por dia. Se multiplica al final para transformar a radianes por minuto 14 n0<- asteRisk:::semiMajorAxisToMeanMotion(koe$semiMajorAxis)*((2*pi)/( 1440)) 15 e0<- koe$eccentricity 16 i0<- koe$inclination 17 M0<- koe$meanAnomaly 18 omega0<- koe$argumentPerigee 19 OMEGA0<- koe$longitudeAscendingNode 20 return(c(n0, e0, i0, M0, omega0, OMEGA0)) 21 } Extracto de c´odigo 7.3: Obtenci´on de par´ametros orbitales a partir de posici´on y velocidad en GCRF 7.2.3. Calcular el error medio de predicci´on Puesto que la idea principal del proyecto es mejorar la precisi´on de SGDP4 frente a HPOP, resulta natural calcular el error cometido con SGDP4 en la predicci´on. La primera implementaci´on fue una funci´on que recib´ıa dos conjuntos de posiciones, pareadas seg´un la ´epoca, es decir, el instante de tiempo, y calculaba la distancia media entre las posiciones. La utilidad inicial era utilizarla con una lista de predicciones de HPOP y otra de SGDP4, aunque finalmente no se utiliza´o para HPOP sino que sirvi´o de ayuda en la interpretaci´on de los resultados de los modelos, para ver el como variaba el error antes y despu´es de aplicar la correcci´on de la red. 1calcula_error_medio <- function(posiciones_reales, posiciones_predichas) { 2 3error_prediccion <- 0 4num_posiciones <- nrow(posiciones_reales) 46 7.2.6. Creaci´on del primer dataset La funci´on que se implement´o devolv´ıa un dataset que incluye las posiciones y las velocidades de las predicciones realizadas por HPOP y SGDP4 a partir de una fecha inicial en una serie de tiempos. Sin embargo, debido al coste computacional que tiene ejecutar HPOP, se vi´o que no resultaba viable generar los datos de esta forma y se decidi´o utilizar los ficheros de datos RINEX que, adem´as, contienen las posiciones y velocidades exactas, salvo peque˜nos errores de observaci´on. 7.2.7. Obtenci´on de datos RINEX Para la obtenci´on de los ficheros de datos RINEX se realiz´o web scrapping, con la librer´ıa BeautifulSoup en Python, ya que hab´ıa multitud de ficheros y descargarlos a mano no era factible. A continuaci´on se muestra la estructura de directorios del repositorio del Instituto de Estad´ıstica y Cartograf´ıa de la Junta de Andaluc´ıa. Dentro de cada fichero ZIP se encuentran cuatro archivos: un archivo html especificando las caracter´ısticas de la antena receptora de las se˜nales, el archivo de datos meteorol´ogicos y los archivos de observaciones GPS y GLONASS. Figura 7.1: Estructura de directorios del repositorio de RINEX versi´on 2 de la Junta de Andaluc´ıa 1import requests 2from bs4import BeautifulSoup 53 3from datetime import date, timedelta 4from zipfile import ZipFile 5from tqdm import tqdm 6from pathlib import Path 7from os import rename, remove 8 9main_path = Path(__file__).parent.resolve() 10 11 year = 2022 12 start_date = date(year, 1,1) 13 end_date = date(year+1,1,1) 14 15 BASE_URL = "https://www.juntadeandalucia.es/ institutodeestadisticaycartografia/rinex/CRDB" 16 17 def find_zip_link_in_url(url:str): 18 reqs = requests.get(url) 19 soup = BeautifulSoup(reqs.text, ’html.parser’) 20 21 for link in soup.find_all(’a’): 22 u = link.get(’href’) 23 if u[-4:] == ".zip": 24 return u 25 26 def daterange(start_date, end_date): 27 for nin range(int((end_date - start_date).days)): 28 yield start_date + timedelta(n) 29 30 for date in tqdm(daterange(start_date, end_date)): 31 day = str(date.day) 32 month = str(date.month) 33 34 day = f"0{day}" if len(day) == 1else day 35 month = f"0{month}" if len(month) == 1else month 36 37 url = f"{BASE_URL}/{year}/{month}/{day}/24H_30seg" 38 39 t = find_zip_link_in_url(url) 40 41 if(t is None): 42 with open(f"{main_path}/exceptions.txt","a") as excep_file: 43 excep_file.write(f"problem with date: {day}/{month}/{year}, url is: {url} \n") 44 continue 45 46 r = requests.get(t, allow_redirects=True) 47 zip_dest = f’./temp/{day}-{month}-{year}.zip’ 54 48 open(zip_dest, ’wb’).write(r.content) 49 50 with ZipFile(zip_dest, ’r’) as zip: 51 filelist = zip.namelist() 52 for fin filelist: 53 if f[-1:] == "g": 54 zip.extract(f, f"{main_path}/downloads") 55 rename(f"{main_path}/downloads/{f}", f"{main_path}/ downloads/{day}-{month}-{year}.g") 56 57 elif f[-1:] == "n": 58 zip.extract(f, f"{main_path}/downloads/") 59 rename(f"{main_path}/downloads/{f}", f"{main_path}/ downloads/{day}-{month}-{year}.n") 60 61 remove(zip_dest) Extracto de c´odigo 7.10: Obtenci´on ed ficheros RINEX de la Junta de Andaluc´ıa Cabe mencionar que las cabeceras de archivos RINEX contienen datos importantes para su interpretaci´on, tales como correcciones temporales que dependen de m´ultiples factores y cambian entre archivos, por lo que no se pueden unificar varios archivos en uno y deben ser almacenados y tratados por separado. Con el c´odigo anterior, cambiando el valor de year, obtenemos los archivos de navegaci´on RINEX de sat´elites GPS y GLONASS de la estaci´on del Instituto de Estad´ıstica y Cartograf´ıa de la Junta de Andaluc´ıa ubicada en de C´ordoba, de los a˜nos 2019 a 2022. 7.2.8. An´alisis y procesado de los datos Una vez descargados todos los archivos RINEX, debemos proceder a su an´alisis, con el fin de determinar cu´al es el sat´elite del que m´as datos se dispone y utilizarlo para entrenar la red neuronal. El an´alisis se realiza leyendo todos los archivos y contando cu´antos mensajes RINEX de cada sat´elite hay. Debido a que las funciones de lectura del paquete devuelven valores distintos seg´un sea GPS o GLONASS tenemos que implementar funciones distintas para cada constelaci´on. Primero necesitamos una funci´on auxiliar que devuelve un vector con los c´odigos de los sat´elites que aparecen en un determinado archivo, repetidos seg´un el n´umero de mensajes en los que se encuentre. La funci´on que cuente el n´umero de mensajes por sat´elite har´a uso de lo anterior. 1get_satellites_GPS_archivo <- function(rinexData){ 2count <- 0 3prncodes <- c() 4for (i in 1:length(rinexData$messages)) { 55 5prncodes <- c(prncodes, rinexData$messages[[i]]$satellitePRNCode) 6} 7return(prncodes) 8} 9 10 get_num_satellites_GPS_ficheros <- function(ruta_base, ficheros){ 11 satelitesGPS <- c() 12 for (f in ficheros){ 13 print(f) 14 faux <- paste0(ruta_base, f) 15 if (file.exists(faux)) { 16 aux <- readGPSNavigationRINEX(faux) 17 satelitesGPS <- c(satelitesGPS, get_satellites_GPS_archivo(aux)) 18 }else { 19 print("NO EXISTE") 20 } 21 } 22 return(table(satelitesGPS)) 23 } Extracto de c´odigo 7.11: C´alculo de n´umero de mensajes RINEX por sat´elites GPS 1get_satellites_GLONASS_archivo <- function(rinexData){ 2satelliteNumbers <- c() 3for (i in 1:length(rinexData$messages)) { 4satelliteNumbers <- c(satelliteNumbers, rinexData$messages[[i]]$ satelliteNumber) 5} 6return(satelliteNumber) 7} 8 9get_num_satellites_GLONASS_ficheros <- function(ruta_base, ficheros){ 10 satelitesGLONASS <- c() 11 12 for (f in ficheros){ 13 print(f) 14 faux <- paste0(ruta_base, f) 15 if (file.exists(faux)) { 16 aux <- readGLONASSNavigationRINEX(faux) 17 satelitesGLONASS <- c(satelitesGLONASS, get_satellites_GLONASS_ archivo(aux)) 18 }else { 19 print("NO EXISTE") 20 } 21 } 22 return(table(satelitesGLONASS)) 23 } Extracto de c´odigo 7.12: C´alculo de n´umero de mensajes RINEX por sat´elites 56 GLONASS Despu´es de ver que el sat´elite con m´as datos era el Kosmos 2514, de la constelaci´on GLONASS, con c´odigo identificativo 17, se procedi´o a leer de nuevo todos los archivos y almacenar la posici´on, velocidad, aceleraci´on y ´epoca de los mensajes de navegaci´on RINEX. Para ello se implement´o la siguiente funci´on: 1get_datos_GLONASS <- function(ruta_base, ficheros, codigo_satelite){ 2dataset <- data.frame("X" =numeric(), 3"Y" =numeric(), 4"Z" =numeric(), 5"dx" =numeric(), 6"dy" =numeric(), 7"dz" =numeric(), 8"ephemerisUTCTime" =character()) 9for (f in ficheros){ 10 print(f) 11 faux <- paste0(ruta_base, f) 12 if (file.exists(faux)) { 13 rinexData <- readGLONASSNavigationRINEX(faux) 14 for (i in 1:length(rinexData$messages)) { 15 satelliteData <- rinexData$messages[[i]] 16 if (satelliteData$satelliteNumber==codigo_satelite){ 17 position_ITRF <- c(satelliteData$positionX, satelliteData$ positionY, satelliteData$positionZ) 18 velocity_ITRF <- c(satelliteData$velocityX, satelliteData$ velocityY, satelliteData$velocityZ) 19 epoch <- as.character.Date(as.POSIXct(satelliteData$ ephemerisUTCTime, tz="UTC")) 20 coordinates_GCRF <- ITRFtoGCRF(position_ITRF, velocity_ITRF, epoch) 21 new_row <- data.frame("X" = coordinates_GCRF$position[1]*1000, 22 "Y" = coordinates_GCRF$position[2]*1000, 23 "Z" = coordinates_GCRF$position[3]*1000, 24 "dX" = coordinates_GCRF$velocity[1]*1000, 25 "dY" = coordinates_GCRF$velocity[2]*1000, 26 "dZ" = coordinates_GCRF$velocity[3]*1000, 27 "ephemerisUTCTime" =as.character.Date(as .POSIXct(satelliteData$ephemerisUTCTime, tz="UTC"))) 28 dataset <- rbind(dataset, new_row) 29 } 30 } 31 }else { 32 print("NO EXISTE") 33 } 34 } 35 return(dataset) 36 } 57 Extracto de c´odigo 7.13: Obtenci´on de datos de sat´elite GLONASS por c´odigo identificativo 7.2.9. Creaci´on del dataset con las predicciones Debido a la optimizaci´on de los par´ametros orbitales, se implementaron varias funciones auxiliares, que posteriormente fueron utilizadas en la funci´on principal que genera el conjunto de datos. Tambi´en se hizo uso de las librer´ıas dplyr ytidyr. Todos los datasets parciales que se iban construyendo, as´ı como el definitivo, se fueron guardando en archivos de texto, para una lectura f´acil y r´apida con la funci´on read.table. Debido a la gran cantidad de columnas del dataset, se implement´o una funci´on para crear un dataframe vac´ıo, con los campos que iba a tener: 1create_clean_df <- function(){ 2return( 3data.frame("UTCTime_base" =character(), 4"unixTime_base" =numeric(), 5"X_base_real" =numeric(), 6"Y_base_real" =numeric(), 7"Z_base_real" =numeric(), 8"dX_base_real" =numeric(), 9"dY_base_real" =numeric(), 10 "dZ_base_real" =numeric(), 11 "time_to_predict" =numeric(), 12 "UTCTime_predict" =character(), 13 "unixTime_prediict" =numeric(), 14 "X_predict_real" =numeric(), 15 "Y_predict_real" =numeric(), 16 "Z_predict_real" =numeric(), 17 "dX_predict_real" =numeric(), 18 "dY_predict_real" =numeric(), 19 "dZ_predict_real" =numeric(), 20 "X_predict_sgdp4_optim" =numeric(), 21 "Y_predict_sgdp4_optim" =numeric(), 22 "Z_predict_sgdp4_optim" =numeric(), 23 "dX_predict_sgdp4_optim" =numeric(), 24 "dY_predict_sgdp4_optim" =numeric(), 25 "dZ_predict_sgdp4_optim" =numeric() 26 ) 27 ) 28 } Extracto de c´odigo 7.14: Creaci´on del dataframe base 58 Otra funci´on implementada es para buscar la fecha m´as cercana a una dada en los archivos RINEX, dentro de un margen de error. Esto es necesario porque en dichos archivos s´olo disponemos de un n´umero finito de valores, aunque pueda ser muy grande. As´ı, si queremos realizar una predicci´on en una ´epoca concreta, lo m´as probable es que no dispongamos de datos de dicha ´epoca. Cuando esto ocurra, deberemos obtener los de la ´epoca m´as cercana disponible, siempre que esa ´epoca no sea muy lejana a la primera. La fecha que se le pasa como par´ametro tiene que ser de tipo car´acter, para evitar problemas con difftime. 1busca_fecha <- function(dataset, fecha, margen_error_fecha){ 2diferencias_tiempos <- abs(as.POSIXct(dataset$ephemerisUTCTime, tz=’UTC ’)-as.POSIXct(fecha, tz=’UTC’)) 3fecha_siguiente <- dataset$ephemerisUTCTime[which.min(diferencias_ tiempos)] 4dif_fecha_siguiente <- difftime(fecha_siguiente, fecha, units="secs") 5if (dif_fecha_siguiente < margen_error_fecha) { 6return(fecha_siguiente) 7}else { 8return("MARGEN SUPERADO") 9} 10 } Extracto de c´odigo 7.15: B´usqueda de fecha m´as cercana en un dataset Tambi´en se desarroll´o una funci´on para crear un dataframe con las predicciones realizadas con SGDP4 a partir de una sola fecha en una serie de tiempos indicados, para utilizarla despu´es en la funci´on para crear el conjunto de datos completo. Tambi´en recibe como par´ametros un dataset base con las datos exactos para poder optimizar los par´ametros, un margen de error para la b´usqueda de la fecha m´as cercana al tiempo de predicci´on, el n´umero de puntos con los que aproximar, adem´as del de la fecha base, y el coeficiente de rozamiento Bstar. Se ha tenido que a˜nadir un manejador de excepciones, ya que en algunas ocasiones se llegaba a valores de los par´ametros que, aunque v´alidos, no eran admisibles por la implementaci´on de la funci´on sgdp4 de asteRisk para realizar predicciones a largo plazo con ellos. 1get_df_predictions_from_one_point_sgdp4<- function(dataset, fecha_base, elementos_orbitales, lista_tiempos, margen_error, puntos_aprox, Bstar) { 2 3## fecha_base es una fecha del dataset 4 5#Inicializamos el dataframe a devolver 6res_df <- create_clean_df() 7 8#Obtenemos los datos para dicha fecha 9data_base <- dataset[dataset$ephemerisUTCTime==fecha_base,][1,] 10 #Nos quedamos con la primera fila en caso de repetidos 11 12 fechas_a_predecir <- c() 59 13 for (tin lista_tiempos){ 14 #Calculamos la fecha en la que vamos a generar la prediccion 15 fecha_a_buscar <- as.character.Date(as.POSIXct(fecha_base, tz=’UTC’)+ t) 16 margen_error_fecha <- (1+margen_error)*t 17 18 #Si dicha fecha no se encuentra en el dataset cogemos la primera siguiente, siempre que no se distancie mucho 19 siguiente_fecha <- busca_fecha(dataset, fecha_a_buscar, margen_error_ fecha) 20 21 if (siguiente_fecha!="MARGEN SUPERADO") { 22 fechas_a_predecir <- c(fechas_a_predecir, siguiente_fecha) 23 } 24 } 25 26 fechas_a_predecir <- fechas_a_predecir[!duplicated(fechas_a_predecir)] 27 28 for (fecha_prediccion in fechas_a_predecir) { 29 tryCatch({ 30 point_TEME <- sgdp4(n0=elementos_orbitales[1], e0=elementos_ orbitales[2], i0=elementos_orbitales[3], 31 M0=elementos_orbitales[4], omega0=elementos_ orbitales[5], 32 OMEGA0=elementos_orbitales[6], Bstar=elementos_ orbitales[7], 33 initialDateTime=fecha_base, targetTime=fecha_ prediccion) 34 point_GCRF <- TEMEtoGCRF(point_TEME$position*1000, point_TEME$ velocity*1000, fecha_prediccion) 35 }, 36 error=function(e){ 37 point_GCRF <- list(NA, NA) 38 names(point_GCRF) <- c("position","velocity") 39 }, 40 warning=function(w){ 41 point_GCRF <- list(NA, NA) 42 names(point_GCRF) <- c("position","velocity") 43 }) 44 45 data_predicted_real <- dataset[dataset$ephemerisUTCTime==fecha_ prediccion, ] 46 47 new_row <- data.frame("ephemerisUTCTime_base" = fecha_base, 48 "unixTime_base" =as.numeric(as.POSIXct(fecha_ base, tz="UTC")), 49 "X_base_real" =data_base$X, 60 50 "Y_base_real" =data_base$Y, 51 "Z_base_real" =data_base$Z, 52 "dX_base_real" =data_base$dX, 53 "dY_base_real" =data_base$dY, 54 "dZ_base_real" =data_base$dZ, 55 "time_to_predicts_mins" =abs(as.POSIXct(fecha_ base, tz=’UTC’)-as.POSIXct(fecha_prediccion, tz=’UTC’)), 56 "ephemerisUTCTime_predict" = fecha_prediccion, 57 "unixTime_predict" =as.numeric(as.POSIXct( fecha_prediccion, tz="UTC")), 58 "X_predict_real" =data_predicted_real$X, 59 "Y_predict_real" =data_predicted_real$Y, 60 "Z_predict_real" =data_predicted_real$Z, 61 "dX_predict_real" =data_predicted_real$dX, 62 "dY_predict_real" =data_predicted_real$dY, 63 "dZ_predict_real" =data_predicted_real$dZ, 64 "X_predict_sgdp4_optim" = point_GCRF$position[1 ], 65 "Y_predict_sgdp4_optim" = point_GCRF$position[2 ], 66 "Z_predict_sgdp4_optim" = point_GCRF$position[3 ], 67 "dX_predict_sgdp4_optim" = point_GCRF$velocity[ 1], 68 "dY_predict_sgdp4_optim" = point_GCRF$velocity[ 2], 69 "dZ_predict_sgdp4_optim" = point_GCRF$velocity[ 3] 70 ) 71 res_df <- rbind(res_df,new_row) 72 } 73 return(res_df) 74 } Extracto de c´odigo 7.16: Creaci´on de dataset a partir de una ´unica fecha Con las funciones anteriores ya podemos implementar la funci´on que genera el conjunto de datos a partir de los par´ametros de entrada iniciales, que son el dataset de las posiciones (obtenido de los ficheros RINEX), un n´umero de puntos que indica cu´antas fechas se tomar´an como base para realizar las predicciones en los periodos indicados en la lista de tiempos, un margen de error para la b´usqueda de las fechas, el n´umero de puntos empleado en la aproximaci´on y el coeficiente de rozamiento, que por defecto estar´a a 0 (ya que no lo conocemos en la mayor´ıa de los casos). 1get_df_predictions_sgdp4<- function(dataset, num_puntos, lista_tiempos, margen_error, puntos_aprox, Bstar=0){ 2 3puntos_aleatorios <- sample(1:nrow(dataset), num_puntos, replace=FALSE ) 61 4fechas_base <- dataset[puntos_aleatorios,]$ephemerisUTCTime 5 6res_df <- create_clean_df() 7 8for (i in 1:length(fechas_base)) { 9print(paste("Punto",i,"de", num_puntos)) 10 fecha <- fechas_base[i] 11 print(fecha) 12 elementos_orbitales <- aprox_params_orbit(dataset=dataset, 13 fecha=fecha, 14 num_puntos=puntos_aprox, 15 Bstar=0) 16 res_parcial <- get_df_predictions_from_one_point_sgdp4(dataset = dataset, fecha_base = fecha, elementos_orbitales = elementos_orbitales , lista_tiempos = lista_tiempos, margen_error = margen_error, puntos_ aprox = puntos_aprox, Bstar = Bstar) 17 res_df <- rbind(res_df, res_parcial) 18 } 19 return(res_df) 20 } Extracto de c´odigo 7.17: Creaci´on de dataset El valor por defecto del par´ametro Bstar se ha establecido a 0, ya que en la mayor´ıa de los casos no vamos a conocerlo y as´ı, al optimizar se acercar´a al valor real. Es posible que en el conjunto de datos generado con la funci´on anterior aparezca alg´un NA, por lo explicado anteriormente. Se decidi´o guardar esas observaciones faltantes por si en alg´un momento se quiere realizar un an´alisis acerca de cu´ando falla la funci´on del paquete. Por ejemplo, podr´ıa interesar estimar cu´al deber´ıa ser la lista de tiempos a predecir para maximizar el n´umero de datos sin que haya problemas con la implementaci´on del modelo SGDP4. De todos modos, para eliminar las filas con observaciones faltantes, bast´o con utilizar la funci´on na.omit de R. Cabe mencionar que se almacen´o en el dataset un campo time to predict, que indica el tiempo de diferencia con el que se busc´o esa fecha en el conjunto de datos anterior. La finalidad de esta variable es por si hiciera falta depurar algo. 7.2.10. Adici´on de las predicciones con par´ametros sin optimizar Las predicciones con los par´ametros optimizados reducen el error a medio plazo, ya que se ajustan los par´ametros a un fragmento de ´orbita m´as extensa, pero lo aumentan a corto, lo cual puede ser se˜nal de que la ´orbita se modifica demasiado al cambiar los par´ametros. Por tanto, tambi´en se incluyeron en el dataset las posiciones predichas por SGDP4 con los par´ametros sin optimizar, para evaluar el comportamiento de la red neuronal. 62 17 18 summary(model) Extracto de c´odigo 7.29: Creaci´on del modelo y resumen del mismo Para el entrenamiento de los mismos se utilizaron un tama˜no de batch de 128 y el n´umero de epochs ser´a de 100, despu´es de haber realizado pruebas con otras combinaciones de par´ametros y observar que esta era la que mejor funcionaba. La funci´on de p´erdida escogida fue el error cuadr´atico medio, mientras que la m´etrica fue el error absoluto medio. 1#Compilacion del modelo 2model %> % 3compile( 4optimizer = "adam", 5loss = "mse",#"mean_squared_error", 6metrics = "mae" #"mean_absolute_error" 7) 8 9#Entrenamiento del modelo 10 history <- model %> % 11 fit( 12 x=x_train_scaled, 13 y=y_train_scaled, 14 batch_size = 128, 15 epochs = 100, 16 verbose = 2 17 ) Extracto de c´odigo 7.30: Compilaci´on y entrenamiento del modelo 7.2.15. Obtenci´on de los datos de los astros Como se coment´o en la parte de dise˜no, se utiliza la funci´on JPLephemerides del paquete asteRisk. El procedimiento es muy similar al de la generaci´on de los otros conjuntos de datos. Como se puede observar no se incluyeron los datos de la Tierra, puesto que al estar utilizando el sistema de referencia GCRF, las coordenadas de la posici´on y las componentes de la velocidad son todas 0. 1datos <- read.table("traintest.txt", header=TRUE) 2 3datos_astros <- data.frame() 4 5fechas_sin_repetir <- datos$ephemerisUTCTime_base[!duplicated(datos$ ephemerisUTCTime_base)] 6 69 7for (fecha in fechas_sin_repetir) { 8MJD_UTC <- dateTimeToMJD(fecha, timeSystem = "UTC") 9ephemerides <- JPLephemerides(MJD_UTC, timeSystem = "UTC", centralBody= "Earth") 10 11 new_row <- data.frame("UTCTime" = fecha, 12 "unixTime" =as.numeric(as.POSIXct(fecha, tz="UTC ")), 13 "X_Sun" = ephemerides$positionSun[1], 14 "Y_Sun" = ephemerides$positionSun[2], 15 "Z_Sun" = ephemerides$positionSun[3], 16 "dX_Sun" = ephemerides$velocitySun[1], 17 "dY_Sun" = ephemerides$velocitySun[2], 18 "dZ_Sun" = ephemerides$velocitySun[3], 19 "ddX_Sun" = ephemerides$accelerationSun[1], 20 "ddY_Sun" = ephemerides$accelerationSun[2], 21 "ddZ_Sun" = ephemerides$accelerationSun[3], 22 "X_SunSSBarycentric" = ephemerides$ positionSunSSBarycentric[1], 23 "Y_SunSSBarycentric" = ephemerides$ positionSunSSBarycentric[2], 24 "Z_SunSSBarycentric" = ephemerides$ positionSunSSBarycentric[3], 25 "dX_SunSSBarycentric" = ephemerides$ velocitySunSSBarycentric[1], 26 "dY_SunSSBarycentric" = ephemerides$ velocitySunSSBarycentric[2], 27 "dZ_SunSSBarycentric" = ephemerides$ velocitySunSSBarycentric[3], 28 "ddX_SunSSBarycentric" = ephemerides$ accelerationSunSSBarycentric[1], 29 "ddY_SunSSBarycentric" = ephemerides$ accelerationSunSSBarycentric[2], 30 "ddZ_SunSSBarycentric" = ephemerides$ accelerationSunSSBarycentric[3], 31 "X_Mercury" = ephemerides$positionMercury[1], 32 "Y_Mercury" = ephemerides$positionMercury[2], 33 "Z_Mercury" = ephemerides$positionMercury[3], 34 "dX_Mercury" = ephemerides$velocityMercury[1], 35 "dY_Mercury" = ephemerides$velocityMercury[2], 36 "dZ_Mercury" = ephemerides$velocityMercury[3], 37 "ddX_Mercury" = ephemerides$accelerationMercury[1 ], 38 "ddY_Mercury" = ephemerides$accelerationMercury[2 ], 39 "ddZ_Mercury" = ephemerides$accelerationMercury[3 ], 70 40 "X_Venus" = ephemerides$positionVenus[1], 41 "Y_Venus" = ephemerides$positionVenus[2], 42 "Z_Venus" = ephemerides$positionVenus[3], 43 "dX_Venus" = ephemerides$velocityVenus[1], 44 "dY_Venus" = ephemerides$velocityVenus[2], 45 "dZ_Venus" = ephemerides$velocityVenus[3], 46 "ddX_Venus" = ephemerides$accelerationVenus[1], 47 "ddY_Venus" = ephemerides$accelerationVenus[2], 48 "ddZ_Venus" = ephemerides$accelerationVenus[3], 49 "X_Mars" = ephemerides$positionMars[1], 50 "Y_Mars" = ephemerides$positionMars[2], 51 "Z_Mars" = ephemerides$positionMars[3], 52 "dX_Mars" = ephemerides$velocityMars[1], 53 "dY_Mars" = ephemerides$velocityMars[2], 54 "dZ_Mars" = ephemerides$velocityMars[3], 55 "ddX_Mars" = ephemerides$accelerationMars[1], 56 "ddY_Mars" = ephemerides$accelerationMars[2], 57 "ddZ_Mars" = ephemerides$accelerationMars[3], 58 "X_Jupiter" = ephemerides$positionJupiter[1], 59 "Y_Jupiter" = ephemerides$positionJupiter[2], 60 "Z_Jupiter" = ephemerides$positionJupiter[3], 61 "dX_Jupiter" = ephemerides$velocityJupiter[1], 62 "dY_Jupiter" = ephemerides$velocityJupiter[2], 63 "dZ_Jupiter" = ephemerides$velocityJupiter[3], 64 "ddX_Jupiter" = ephemerides$accelerationJupiter[1 ], 65 "ddY_Jupiter" = ephemerides$accelerationJupiter[2 ], 66 "ddZ_Jupiter" = ephemerides$accelerationJupiter[3 ], 67 "X_Saturn" = ephemerides$positionSaturn[1], 68 "Y_Saturn" = ephemerides$positionSaturn[2], 69 "Z_Saturn" = ephemerides$positionSaturn[3], 70 "dX_Saturn" = ephemerides$velocitySaturn[1], 71 "dY_Saturn" = ephemerides$velocitySaturn[2], 72 "dZ_Saturn" = ephemerides$velocitySaturn[3], 73 "ddX_Saturn" = ephemerides$accelerationSaturn[1], 74 "ddY_Saturn" = ephemerides$accelerationSaturn[2], 75 "ddZ_Saturn" = ephemerides$accelerationSaturn[3], 76 "X_Uranus" = ephemerides$positionUranus[1], 77 "Y_Uranus" = ephemerides$positionUranus[2], 78 "Z_Uranus" = ephemerides$positionUranus[3], 79 "dX_Uranus" = ephemerides$velocityUranus[1], 80 "dY_Uranus" = ephemerides$velocityUranus[2], 81 "dZ_Uranus" = ephemerides$velocityUranus[3], 82 "ddX_Uranus" = ephemerides$accelerationUranus[1], 83 "ddY_Uranus" = ephemerides$accelerationUranus[2], 71 84 "ddZ_Uranus" = ephemerides$accelerationUranus[3], 85 "X_Neptune" = ephemerides$positionNeptune[1], 86 "Y_Neptune" = ephemerides$positionNeptune[2], 87 "Z_Neptune" = ephemerides$positionNeptune[3], 88 "dX_Neptune" = ephemerides$velocityNeptune[1], 89 "dY_Neptune" = ephemerides$velocityNeptune[2], 90 "dZ_Neptune" = ephemerides$velocityNeptune[3], 91 "ddX_Neptune" = ephemerides$accelerationNeptune[1 ], 92 "ddY_Neptune" = ephemerides$accelerationNeptune[2 ], 93 "ddZ_Neptune" = ephemerides$accelerationNeptune[3 ], 94 "X_Pluto" = ephemerides$positionPluto[1], 95 "Y_Pluto" = ephemerides$positionPluto[2], 96 "Z_Pluto" = ephemerides$positionPluto[3], 97 "dX_Pluto" = ephemerides$velocityPluto[1], 98 "dY_Pluto" = ephemerides$velocityPluto[2], 99 "dZ_Pluto" = ephemerides$velocityPluto[3], 100 "ddX_Pluto" = ephemerides$accelerationPluto[1], 101 "ddY_Pluto" = ephemerides$accelerationPluto[2], 102 "ddZ_Pluto" = ephemerides$accelerationPluto[3], 103 "X_Moon" = ephemerides$positionMoon[1], 104 "Y_Moon" = ephemerides$positionMoon[2], 105 "Z_Moon" = ephemerides$positionMoon[3], 106 "dX_Moon" = ephemerides$velocityMoon[1], 107 "dY_Moon" = ephemerides$velocityMoon[2], 108 "dZ_Moon" = ephemerides$velocityMoon[3], 109 "ddX_Moon" = ephemerides$accelerationMoon[1], 110 "ddY_Moon" = ephemerides$accelerationMoon[2], 111 "ddZ_Moon" = ephemerides$accelerationMoon[3], 112 "lunar_libration_Phi" = ephemerides$ lunarLibrationAngles[1], 113 "lunar_libration_Theta" = ephemerides$ lunarLibrationAngles[2], 114 "lunar_libration_Psi" = ephemerides$ lunarLibrationAngles[3], 115 "lunar_libration_dPhi" = ephemerides$ lunarLibrationAnglesDerivatives[1], 116 "lunar_libration_dTheta" = ephemerides$ lunarLibrationAnglesDerivatives[2], 117 "lunar_libration_dPsi" = ephemerides$ lunarLibrationAnglesDerivatives[3], 118 "lunar_libration_ddPhi" = ephemerides$ lunarLibrationAnglesSecondDerivatives[1], 119 "lunar_libration_ddTheta" = ephemerides$ lunarLibrationAnglesSecondDerivatives[2], 72 120 "lunar_libration_ddPsi" = ephemerides$ lunarLibrationAnglesSecondDerivatives[3] 121 ) 122 datos_astros <- rbind(datos_astros, new_row) 123 } 124 125 write.table(datos_astros, file ="datosastros.txt",row.names = FALSE) Extracto de c´odigo 7.31: Obtenci´on de los datos de los astros 7.2.16. Adici´on de los datos de los astros al dataset Una vez que tenemos los datos de los astros del Sistema Solar procedemos a incluirlos en el dataset. El c´odigo empleado es el siguiente: 1datos <- read.table("traintest.txt", header=TRUE) 2datos_astros <- read.table("datosastros.txt", header=TRUE) 3 4datos_ampliados <- data.frame() 5 6for (i in 1:nrow(datos)) { 7# Un contador para tener referencia de en que punto del proceso nos encontramos 8if (i % %1000==0) { 9print(i) 10 } 11 datos_satelite <- datos[i,] 12 datos_astros_satelite <- datos_astros[(datos_astros$UTCTime==datos_ satelite$ephemerisUTCTime_base), -1] 13 new_row <- cbind(datos_satelite, datos_astros_satelite) 14 datos_ampliados <- rbind(datos_ampliados, new_row) 15 } 16 17 write.table(datos_ampliados, file ="traintestconastros.txt",row.names = FALSE) Extracto de c´odigo 7.32: A˜nadir datos de los astros al dataset 7.2.17. Creaci´on de modelos con datos de los astros La ´unica diferencia con la creaci´on de los modelos sin los datos de los astros est´a en la selecci´on de campos del dataframe. A continuaci´on se muestra un ejemplo de selecci´on de campos. 1library(dplyr) 2datos <- select(datos, unixTime_base, 3X_base_real, 73 4Y_base_real, 5Z_base_real, 6dX_base_real, 7dY_base_real, 8dZ_base_real, 9X_Sun, 10 Y_Sun, 11 Z_Sun, 12 dX_Sun, 13 dY_Sun, 14 dZ_Sun, 15 ddX_Sun, 16 ddY_Sun, 17 ddZ_Sun, 18 X_Venus, 19 Y_Venus, 20 Z_Venus, 21 dX_Venus, 22 dY_Venus, 23 dZ_Venus, 24 ddX_Venus, 25 ddY_Venus, 26 ddZ_Venus, 27 X_Moon, 28 Y_Moon, 29 Z_Moon, 30 dX_Moon, 31 dY_Moon, 32 dZ_Moon, 33 ddX_Moon, 34 ddY_Moon, 35 ddZ_Moon, 36 lunar_libration_Phi, 37 lunar_libration_Theta, 38 lunar_libration_Psi, 39 lunar_libration_dPhi, 40 lunar_libration_dTheta, 41 lunar_libration_dPsi, 42 lunar_libration_ddPhi, 43 lunar_libration_ddTheta, 44 lunar_libration_ddPsi, 45 X_Mars, 46 Y_Mars, 47 Z_Mars, 48 dX_Mars, 49 dY_Mars, 50 dZ_Mars, 74 51 ddX_Mars, 52 ddY_Mars, 53 ddZ_Mars, 54 X_Jupiter, 55 Y_Jupiter, 56 Z_Jupiter, 57 dX_Jupiter, 58 dY_Jupiter, 59 dZ_Jupiter, 60 ddX_Jupiter, 61 ddY_Jupiter, 62 ddZ_Jupiter, 63 unixTime_predict, 64 X_predict_sgdp4, 65 Y_predict_sgdp4, 66 Z_predict_sgdp4, 67 dX_predict_sgdp4, 68 dY_predict_sgdp4, 69 dZ_predict_sgdp4, 70 error_X_predict, 71 error_Y_predict, 72 error_Z_predict 73 ) Extracto de c´odigo 7.33: Elecci´on de campos del dataset con los astros 7.2.18. Evaluaci´on de los modelos Para la evaluaci´on de los modelos creados se evaluaron las m´etricas sobre el conjunto de test. Como dichas m´etricas proven´ıan de datos normalizados no eran f´aciles de interpretar, por lo que se procedi´o a la desnormalizaci´on y a calcular el error medio de la predicci´on antes y despu´es de aplicar la correcci´on de la red neuronal. 1result <- model %> % evaluate(x_test, y_test) 2result["mae"] Extracto de c´odigo 7.34: Evaluaci´on de m´etricas sobre el conjunto de test 1predictions_normalized <- model %> % predict(x_test_scaled) Extracto de c´odigo 7.35: Predicciones con la red Recordemos que la desnormalizaci´on tiene que realizarse con los par´ametros obtenidos del conjunto de entrenamiento. 1predictions <- matrix(nrow=nrow(predictions_normalized), ncol=ncol( predictions_normalized)) 2for (i in 1:ncol(predictions_normalized)) { 75 3predictions[,i] <- predictions_normalized[,i]*iqr_y[i] + mediana_y[i] 4} Extracto de c´odigo 7.36: Desnormalizaci´on de las predicciones Para la correcci´on de las predicciones de SGDP4 basta sumarle la predicci´on realizada por la red, desnormalizada. 1predicciones_sin_corregir <- select(df_test, X_predict_sgdp4_optim, 2Y_predict_sgdp4_optim, 3Z_predict_sgdp4_optim) 4colnames(predicciones_sin_corregir) <- c("X","Y","Z") 5 6 7predicciones_corregidas <- data.frame("X" =numeric(), 8"Y" =numeric(), 9"Z" =numeric()) 10 11 for (i in 1:nrow(predicciones_sin_corregir)) { 12 new_row <- data.frame("X" = predicciones_sin_corregir[i,]$X + error_ predictions_df[i,]$X, 13 "Y" = predicciones_sin_corregir[i,]$Y + error_ predictions_df[i,]$Y, 14 "Z" = predicciones_sin_corregir[i,]$Z + error_ predictions_df[i,]$Z) 15 predicciones_corregidas <- rbind(predicciones_corregidas, new_row) 16 } Extracto de c´odigo 7.37: Correcci´on de predicci´on SGDP4 1real_test <- select(df_test, X_predict_real, 2Y_predict_real, 3Z_predict_real) 4colnames(real_test) <- c("X","Y","Z") Extracto de c´odigo 7.38: Obtenci´on de posiciones reales en el conjunto de test De este modo ya podemos hacer uso de la funci´on que recib´ıa dos conjuntos con posiciones y calculaba el error medio entre las pareadas. 7.3. Conclusiones En este cap´ıtulo se han expuesto los detalles de la implementaci´on, documentando las funciones empleadas para todas las tareas, de modo que los avances realizados puedan ser utilizados en el futuro. De este modo, para la optimizaci´on de par´ametros y la creaci´on de nuevos conjuntos de datos de sat´elites se podr´a reutilizar el c´odigo, ahorrando tiempo y costes. 76 8. Pruebas y experimentos En este cap´ıtulo se tratar´an los modelos entrenados, explicando qu´e par´ametros reciben y evaluando la validez y eficacia de su uso. Cuando se mencione que un modelo recibe una posici´on, un movimiento o una velocidad, nos estaremos refiriendo a que recibe como entrada las 3 coordenadas o componentes de dicha magnitud (X, Y y Z), an´alogamente con la ostras magnitudes como la libraci´on lunar. Las fechas se considerar´an en formato UNIX. Todos los modelos tuvieron como par´ametros de entrada la fecha base de la predicci´on, la posici´on del sat´elite en dicha fecha y la fecha en la que realizar la predicci´on. As´ı, nos referimos a ellos como par´ametros base. Separaremos los modelos en dos grandes familias, aquellos que no reciben los datos de los astros y aquellos que s´ı. Dentro de los que no reciben datos de los astros, primero consideramos los modelos que reciben la posici´on predicha y luego pasaremos a aquellos que reciben el movimiento. A su vez, dentro de ellos separaremos los que involucran a las predicciones con los par´ametros optimizados de las que se han realizado sin optimizarlos previamente. Todos los modelos ir´an acompa˜nados de una gr´afica que muestra c´omo ha variado el error absoluto medio (MAE) a lo largo de su entrenamiento. No se muestra la evoluci´on de la p´erdida (loss) durante el entrenamiento porque la funci´on elegida fue el error cuadr´atico medio (MSE), que muestra un comportamiento muy similar al MAE y su interpretaci´on menos intuitiva. Previamente a la evaluaci´on de los modelos comentaremos los errores medios de las predicciones antes de aplicar la correcci´on de las redes neuronales. 8.1. Error medio de las predicciones SGDP4 Se calcul´o el error medio de las predicciones realizadas con SGDP4 con respecto a la posici´on real. Para dicho c´alculo se tuvo en cuenta los datos del conjunto de prueba, de modo que posteriormente se pudiera comparar con el de las predicciones tras corregir el error con las redes neuronales. Los resultados se recogen en la siguiente tabla: Predicci´on SGDP4 Error medio (metros) Sin optimizar elementos orbitales 578619 Con los elementos orbitales optimizados 25022408 Cuadro 8.1: Error medio de las predicciones antes de corregirlas con los resultados de la red Se observa que el error tras optimizar los par´ametros es superior al otro, seguramente es debido a que la optimizaci´on se realiza en el entorno del punto base de la 77 predicci´on, de modo que la ´orbita que mejor se ajuste a ellos no sea la que mejor lo haga a la trayectoria real, ocasionando errores muy grandes en predicciones a largo plazo. 8.2. Evaluaci´on de los modelos sin datos de los astros En primer lugar se analizan los modelos que reciben las predicciones con los par´ametros sin optimizar. 8.2.1. Par´ametros orbitales sin optimizar Modelos que reciben la posici´on predicha El primer modelo que recibe los par´ametros base y la predicci´on realizada con SGDP4 sin optimizar los par´ametros orbitales previamente. Sobre el conjunto de entrenamiento obtuvo un MAE de 13.6329, mientras que sobre el de test fue 14.39973. El error medio de la predicci´on de posiciones despu´es de la correcci´on fue de 647040.5, que empeora la predicci´on. Figura 8.1: Evoluci´on del MAE recibiendo los par´ametros base y la posici´on predicha con los par´ametros orbitales sin optimizar Posteriormente se a˜nadi´o la velocidad en el instante inicial como par´ametro de entrada, obteniendo as´ı un modelo con un MAE de 12.0608 sobre entrenamiento y 12.98055 sobre test. El error medio de la predicci´on de posiciones con la correcci´on fue 78 Figura 8.10: Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitales sin optimizar El modelo que recibe el movimiento predicho, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho sin optimizar obtuvo un MAE de 7.8431 sobre el conjunto de entrenamiento y de 9.702709 sobre el de test. El error medio de la predicci´on tras aplicar la correcci´on fue de 433590.8 metros. Figura 8.11: Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho con los par´ametros orbitales sin optimizar 85 Agrupando estos resultados: Entradas extra del modelo MAE Error medio de la predicci´on corregida Movimiento predicho sin optimizar par´ametros 12.11 540895.5Posiciones de todos los astros del Sistema Solar Libraci´on lunar Movimiento predicho sin optimizar par´ametros 11.55 519358.5 Velocidad inicial Posiciones de todos los astros del Sistema Solar Libraci´on lunar Movimiento predicho sin optimizar par´ametros 9.7 433590.8Velocidad inicial Posiciones del Sol, la Luna y J´upiter Cuadro 8.4: Comparativa de modelos para predicciones con par´ametros sin optimizar y datos de los astros 8.3.2. Par´ametros orbitales optimizados El modelo que recibe el movimiento predicho, las posiciones de todos los astros, la libraci´on lunar y el movimiento predicho optimizado obtuvo un MAE de 0.3509 en entrenamiento y de 0.4538257 en test. El error medio de la predicci´on tras la correcci´on fue de 19314206 metros. Figura 8.12: Evoluci´on del MAE recibiendo los par´ametros base, los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitales optimizados El modelo que recibe el movimiento predicho, la velocidad inicial, las posiciones 86 de todos los astros, la libraci´on lunar y el movimiento predicho optimizado obtuvo un MAE de 0.3532 en entrenamiento y de 0.4563037 en test. El error medio de la predicci´on tras la correcci´on fue de 19404345 metros. Figura 8.13: Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial los datos de todos los astros del Sistema Solar y el movimiento predicho con los par´ametros orbitales optimizados El modelo que recibe el movimiento predicho, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho con los par´ametros optimizados obtuvo un MAE en entrenamiento de 0.3395 y en test de 0.4589353. El error medio tras la predicci´on fue de 19520238 metros. 87 Figura 8.14: Evoluci´on del MAE recibiendo los par´ametros base, la velocidad inicial, las posiciones del Sol, Venus, la Luna, Marte, J´upiter y el movimiento predicho con los par´ametros orbitales optimizados Agrupando estos resultados: Entradas extra del modelo MAE Error medio de la predicci´on corregida Movimiento predicho con par´ametros optim. 0.45 19314206Posiciones de todos los astros del Sistema Solar Libraci´on lunar Movimiento predicho con par´ametros optim. 0.46 19404345 Velocidad inicial Posiciones de todos los astros del Sistema Solar Libraci´on lunar Movimiento predicho con par´ametros optim. 0.46 19520238Velocidad inicial Posiciones del Sol, la Luna y J´upiter Cuadro 8.5: Comparativa de modelos para predicciones con par´ametros optimizados y datos de los astros 8.4. Conclusiones Tras entrenar diversos modelos y analizar el error metido cometido en la predicci´on tras aplicar la correcci´on de los mismos concluimos que es viable el uso de redes neuronales para mejorar la predicci´on del modelo SGDP4. 88 Extradas extra del modelo Par´ametros optimizados MAE Error medio de la predicci´on corregida ∆Error( %) Posici´on predicha No 14.4 647040.5 11.82 Posici´on predicha No 12.98 583126.3 0.78 Velocidad inicial Movimiento predicho No 12.29 559453 -3.31 Movimiento predicho No 13.93 631188.6 9.09 Velocidad inicial Posici´on predicha S´ı 0.41 17331023 -30.74 Posici´on predicha S´ı 0.42 17869264 -28.59 Velocidad inicial Movimiento predicho S´ı 0.39 16831337 -32.73 Movimiento predicho S´ı 0.39 16814065 -32.8 Velocidad inicial Movimiento predicho No 12.11 540895.5 -6.52Posiciones de los astros Libraci´on lunar Movimiento predicho No 11.55 519358.5 -10.24 Velocidad inicial Posiciones de los astros Libraci´on lunar Movimiento predicho No 9.7 433590.8 -25.06 Velocidad inicial Posici´on del Sol Posici´on de la Luna Posici´on de J´upiter Movimiento predicho S´ı 0.45 19314206 -22.81Posiciones de los astros Libraci´on lunar Movimiento predicho S´ı 0.46 19404345 -22.45 Velocidad inicial Posiciones de los astros Libraci´on lunar Movimiento predicho S´ı 0.46 19520238 -21.99 Velocidad inicial Posici´on del Sol Posici´on de la Luna Posici´on de J´upiter Cuadro 8.6: Comparativa de todos los modelos entrenados Destacar el mejor modelo de todos los de las predicciones sin optimizar, que recibe los par´ametros base, la velocidad inicial, las posiciones del Sol, la Luna y J´upiter y el movimiento predicho y que reduce el error medio en un 25 %. 89 En el caso de las predicciones con los par´ametros optimizados, todos los modelos han mejorado la predicci´on, reduciendo considerablemente el error respecto a las predicciones de SGDP4, sin embargo, como estas predicciones eran peores que las otras, la red no ha podido corregir suficientemente el error como para aproximarse a la posici´on real. 90 9. Conclusiones Este Trabajo Fin de Grado se ha centrado en torno a la propagaci´on de trayectorias de sat´elites y la dificultad existente para obtener predicciones precisas, especialmente a medio y largo plazo, con un coste computacional admisible. En este sentido se han desarrollado una serie de algoritmos que permiten optimizar los par´ametros orbitales utilizados por los propagadores SGP4 y SDP4, de modo que sus predicciones se ajusten mejor a la ´orbita en el entorno del punto inicial. Dichos algoritmos fueron tambi´en empleados en la elaboraci´on de un conjunto de datos para el sat´elite Kosmos 2514, que re´une posiciones, velocidades, tanto del sat´elite como de los astros del Sistema Solar, y predicciones sobre la posici´on y velocidad del mismo con el propagador SGDP4, tanto con los par´ametros optimizados como sin optimizar. Cabe mencionar que dichos datos se han unificado para estar en el mismo marco de referencia. Una vez creado el conjunto de datos se utiliz´o para entrenar redes neuronales que predijeran el error cometido con el propagador y as´ı poder corregirlo. Los resultados obtenidos fueron muy prometedores, consiguiendo un modelo que reduc´ıa el error medio en un 25 %. La idea inicial era entrenar tambi´en modelos de redes neuronales recurrentes, que pudieran detectar patrones en los vectores de estado orbitales; sin embargo, debido a la complejidad de la materia, surgieron muchos problemas, como el de la optimizaci´on de par´ametros orbitales que ocup´o m´as tiempo del previsto. A pesar de todo, ha servido para conocer un poco c´omo es el mundo de la investigaci´on, en algunos momentos muy frustrante pero, una vez que se logran los objetivos, muy satisfactorio. Se proponen diversas ideas como trabajo futuro: 1. Entrenar de redes neuronales recurrentes que puedan detectar patrones peri´odicos en los vectores de estado orbital. 2. Integrar con alg´un software de visualizaci´on avanzado. 3. Recopilar m´as datos de fen´omenos astron´omicos que puedan afectar al movimien- to de los sat´elites como, por ejemplo, las tormentas solares. 4. Estudiar en profundidad la optimizaci´on de los par´ametros orbitales para intentar que no ocurra la explosi´on del error de la predicci´on en tiempos lejanos. A t´ıtulo personal, tengo que destacar mi satisfacci´on por haber logrado la elaboraci´on del conjunto de datos, con todas las herramientas matem´aticas y t´ecnicas que han sido necesarias, y por aportar una soluci´on ´util que permita realizar buenas predicciones de trayectorias en un tiempo reducido, aunando todas las competencias adquiridas durante estos cinco a˜nos de Ingenier´ıa Inform´atica y Matem´aticas. 91 10. Bibliograf´ıa [1] Roger R. Bate, Donald D. Mueller, and Jerry E. White. Fundamentals of Astrodynamics. Dover Publications, 1 edition, 1971. [2] Celestrak. Datos actualizados de norad en formato tle. URL https://celestrak. org/NORAD/elements/. [3] Instituto de Estad´ıstica y Cartograf´ıa. Repositorio de ficheros rinex ver- si´on 2 de la junta de andaluc´ıa. URL https://www.juntadeandalucia.es/ institutodeestadisticaycartografia/rinex/. [4] Richard L. Burdenan Douglas J. Faires and Annette M. Burden. An´alisis Num´erico. CENGAGE Learning, 10 edition, 2016. [5] Felix R. Hoods and Ronald L. Roehrich. SPACETRACK REPORT NO. 3. Models for Propagation of NORAD Element Sets. U.S. Department of Defensea, 1988. [6] Richard S. Hujsak. A Restricted Four Body Solution for Resonating Satellites Without Drag. Aerospace Defense Command United Space Air Force, 1979. [7] Zac Yung-Chun Liu, Scott Tarlow, Mohammad Akbar, Quentin Donnellan, and Darek Senkow. Improved orbital propagator integrated with sgp4 and machine learning. In 35th Annual Small Satellite Conference, 2021. URL https://digitalcommons.usu.edu/cgi/viewcontent.cgi?article=5067& context=smallsat. [8] John H. Mathews and Kurtis Fink. M´etodos Num´ericos con MATLAB. Prentice Hall, 3 edition, 2000. [9] Oliver Montenbruck and Eberhard Gill. Satellite Orbits: Models, Methods and Applications. Springer, 1 edition, 2000. [10] Alaa Osama, Mourad Raafat, Ashraf Darwish, Sara Abdelghafar, and Aboul Hassanien. Satellite orbit prediction based on recurrent neural network using two line elements. pages 298–302, 03 2022. doi: 10.1109/ICCI54321.2022.9756063. [11] H. Peng and X. Bai. Limits of Machine Learning Approach on Improving Orbit Prediction Accuracy using Support Vector Machine. In S. Ryan, editor, Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference, page 15, January 2017. [12] Hao Peng and Xiaoli Bai. Artificial neural network–based machine learning approach to improve orbit prediction accuracy. Journal of Spacecraft and Rockets, 55:1–13, 06 2018. doi: 10.2514/1.A34171. [13] Haoli Ren, Xiaolin Chen, Bei Guan, Yongji Wang, Tiantian Liu, and Kongyang Peng. Research on satellite orbit prediction based on neural network algorithm. pages 267–273, 06 2019. ISBN 978-1-4503-7185-8. doi: 10.1145/3341069.3342995. 92 [14] Naciones Unidas. ´ Indice en l´ınea de objetos lanzados al espacio exterior. URL https://www.unoosa.org/oosa/osoindex/search-ng.jspx. [15] David A. Vallado. Fundamentals of Astrodynamics and Applications. Springer, 4 edition, 2013. 93