Analysis and Design of Algorithms for Earth Observation Satellites

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Analysis and Design of Algorithms for Earth Observation Satellites Autor: Manuel Francisco Fernández Melgosa Tutores: Rafael Vázquez Valenzuela, Federico Perea Rojas-Marcos Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Analysis and Design of Algorithms for Earth Observation Satellites Autor: Manuel Francisco Fernández Melgosa Tutores: Rafael Vázquez Valenzuela, Federico Perea Rojas-Marcos Catedrático, Profesor titular Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Trabajo Fin de Máster: Analysis and Design of Algorithms for Earth Observation Satellites Autor: Manuel Francisco Fernández Melgosa Tutores: Rafael Vázquez Valenzuela, Federico Perea Rojas-Marcos El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos El trabajo presentado en este TFM culmina seis años de estudios en ingeniería, un periodo del que me llevo muchas más alegrías que tristezas. Ha sido un camino duro, y en numerosas ocasiones he considerado abandonar al no ver los resultados esperados. Por eso, en primer lugar, quiero expresar mi más profundo agradecimiento a las personas que han hecho posible que llegue hasta aquí: mis padres y mi pareja. Sin su apoyo incondicional, no sé dónde estaría hoy, así que, de corazón, gracias. En relación a este TFM, quiero manifestar toda mi gratitud a los profesores Rafael y Federico, mis tutores, por su orientación constante y sus sinceras recomendaciones sobre mi trabajo. Agradezco también a Deimos Space por la oportunidad laboral que me brindaron mientras realizaba este proyecto. En especial, quiero dar las gracias a Ángel Monge, quien me ha guiado en la empresa y apoyado en la realización de este trabajo, y a Óscar González, por la confianza depositada en mí para este proyecto. Finalmente, quiero dar las gracias a todos mis compañeros de carrera, tanto a aquellos con los que he compartido estos dos últimos años de máster como a los amigos de los primeros años de grado. Gracias a vosotros, las risas después de exámenes suspendidos, presentaciones desastrosas y prácticas fallidas han sido más satisfactorias que los logros académicos. Al final, eso es la ingeniería: cometer errores el 99% del tiempo para que el 1% que cuenta salga perfecto. Gracias de corazón; me llevo una nueva familia para toda la vida. I Resumen Título del trabajo Análisis y Diseño de Algoritmos para Satélites de Observación Terrrestre Resumen del trabajo Los satélites de observación terrestre capturan imágenes y recogen datos del planeta, permitiendo monitorear el medio ambiente, analizar el clima y gestionar desastres, entre otras aplicaciones. La planificación eficiente de su operación es clave debido a la complejidad de coordinar varios satélites en diferentes órbitas y modos de observación. Este trabajo desarrolla algoritmos que permitan planificar estas operaciones de manera automatizada. Partiendo de un modelo teórico, se buscan buenas soluciones para cubrir áreas específicas del planeta. Se proponen algoritmos exactos tipo Branch and Bound, y otros heurísticos y metaheurísticos tipo Greedy y GRASP. La calidad de tales algoritmos se compara con las soluciones obtenidas en un Modelo de Programación Lineal Entera Mixta (MILP por sus siglas en inglés). Finalmente, los experimentos evalúan la eficiencia de los algoritmos en términos de cobertura y tiempo de cálculo, logrando maximizar de manera eficiente la cobertura de las regiones prioritarias. Palabras clave Satélites de observación terrestre, algoritmos de optimización, programación de misiones, Programación Lineal Entera Mixta (MILP), GRASP, Branch and Bound, optimización combinatoria. Conclusiones El objetivo de este trabajo fue desarrollar algoritmos para planificar adquisiciones de satélites, considerando restricciones y preservando las programaciones previas. Se probaron tres enfoques, sin que ninguno fuera el mejor en todas las situaciones probadas. El algoritmo Greedy es el más rápido y ofrece soluciones cercanas a las óptimas. GRASP mejora la cobertura y permite más flexibilidad, aunque con un mayor tiempo de cálculo dependiendo de los parámetros escogidos. Branch and Bound tiene el mayor tiempo de cómputo, pero sigue siendo útil en casos específicos. En última instancia, la elección del mejor algoritmo depende de las características de la misión. III 2Chapter 1. Introduction This is why efficient algorithms that can schedule the pointing of a set of satellites to cover the information of a certain region need to be designed. 1.1 Objectives and motivation The problem addressed in this project is the Earth Observation Satellite Scheduling Problem (EOSSP). The task involves selecting a subset of satellite acquisitions to cover a predefined set of priority regions within a specified time frame, while optimizing a given objective function. Drawing on the reference work [ 5 ], this planning problem can be modelled as a set covering problem, utilizing fundamental mathematical programming techniques. EOS task planing and scheduling is a typical complex combinatorial optimization problem with NP-Hard computational complexity [ 7 ]. Overall, effective task planning and scheduling for Earth Observation Satellites requires a careful balance of mission objectives, technical constraints, and environmental factors to ensure the efficient collection of high-quality data for a range of applications. Advanced algorithms and automation are increasingly being employed to enhance the efficiency and effectiveness of this process. 1.2 Scope of the project The objective of this thesis was proposed by Deimos Space, a Spanish company, to some researchers of the department of Aerospace Engineering and the department of Applied Mathematics II of the University of Seville. Regarding EOSS algorithms, the project aims to explore potential solutions for the scheduling operations problem to gather information about regions of interest that cannot be covered with a single pass of satellites. Deimos Space has previously developed tools for scheduling DEIMOS-2 missions to produce feasible acquisition sequences from a set of user areas of interest [ 10 ]. This topic has been previously studied in various thesis works. In [ 6 ], good results are obtained using heuristic algorithms, such as GRASP without local search, by scheduling acquisitions obtained from satellite constellations with different pointing modes. The work [ 3 ] presented a more complex model of the problem, incorporating constraints such as download links and decomposing the problem into subproblems to explore different solution approaches. In all the referenced works, it is clear that the solution to this problem must be sought combining three different areas: orbital mechanics, computational geometry and algorithms development. Thanks to these previous works, the starting point for this Master’s Thesis is well advanced, particularly in the area of orbital mechanics, which is already highly developed. Significant changes will be made in the block related to computational geometry and algorithms. Thus, along with the requirements of Deimos Space, the objectives of the project are the following: • The study of the systems to calculate the covered zones for the areas defined based on the selected acquisitions. •The calculation of the possible input acquisitions to the selection algorithm. • Based on the calculated possible acquisitions and the areas of interest, the set of acquisitions that maximize or minimize a certain evaluation function will be selected. • The selected acquisitions must meet a set of conditions and available resources. It will be necessary to study the possibility of incorporating these conditions into the acquisition selection algorithm or, alternatively, if these conditions can be checked in a second step by the mission planning system. Intermediate solutions that incorporate a subset of the conditions into the selection algorithm could be feasible. 1.3 Structure of the document 3 • The developed algorithms should allow for the updating of the status of previously selected acquisitions. 1.3 Structure of the document 1. Introduction. Description of the Earth Observation Scheduling Problem (EOSSP) and all the information useful about the realm of the scheduling problems. Scope of the thesis. 2. Problem Setting. Description of all the elements involved in the EOSSP. Formulation of the problem. Optimal solution explained. 3. Data Model. Set of parts conforming the program used to obtain the solution of the problem, with their inputs, outputs and internal structure. 4. Algorithms. Explanation of the different techniques used to solve the EOSSP. 5. Experiments. Set of computational experiments conducted to test the capabilities of the algorithms developed in real situations. 6. Conclusions and Future Work. 2 Problem Setting This chapter introduces the Earth Observation Satellite Scheduling Problem (EOSSP) input data and the formulation of the problem. 2.1 Inputs The initial data of the Earth Observation Satellite Scheduling Problem are as follows: Region A region is a zone of the Earth’s surface and it is defined by its boundary points. Consider as set of regions R={R1,R2,...Rr}, where all regions are pairwise disjoint, that is: Ri∩Rj=/0,i=j,(i,j) = 1,...,r.(2.1) In MatLab, the structure polyshape is used to define any type of polygon [ 9 ], and it is used in this project to define the regions to be covered. This structure has three types of attributes: • Vertices. Points defining the boundary of the region. In case of political zones, this set of points is the frontier. • NumRegions. Number of non-connected parts inside the polyshape structure. According to Equation (2.1), this field can store multiple regions inside it (i.e. multiple islands forming an archipelago), but for the sake of simplifications, each polyshape defining a region has only one element in NumRegions. To work with different regions, each of them is stored in a field of a polyshape array. • NumHoles. Number of holes inside the polyshape structure. Holes are not considered inside the regions. Each region has two parameters: priority and area. The priority pi of a region Ri∈R indicates the importance of the region. If pi>pj , Ri is more important in the planning schedule than Rj (i,j) = 1,...r. In general, the priority of a region is not associated with its area, so a region with greater area could have less priority. The area of each Ri∈Ris denoted by σi. Swaths Consider a satellite with a sensor pointing to the Earth’s surface with different modes. Each mode is an attitude of the satellite, defined by a pitch and roll angles. Swath refers to the area on the Earth’s surface that is imaged or observed by the satellite sensor within a time horizon, defined by its boundary points. Each mode of the satellite generates one swath. In this project, each mode of 5 6Chapter 2. Problem Setting the satellite is considered as a satellite itself, so C satellites with M modes each is equivalent to work with C˙M satellites with one fixed attitude each. Swaths are an input because Deimos Space has tools to propagate satellites orbits and generate the swaths. Acquisitions Consider a swath Wand a region R. An acquisition Ais a set of points inside Wthat covers a part of R, that is, A⊆W/A∩R=/0.(2.2) Acquisitions are stored in polyshape structures, as well as regions. One swath can generate multiple acquisitions depending on the maximum and minimum size of the acquisitions and the size of the region to be covered. To work with the same variables as Deimos Space, each acquisition generated has nine different attributes: •ID Satellite. Satellite used to obtain the acquisition •Attitude. Roll and pitch angles used to generate the acquisition. •Segment. Set of points defining the swath where the acquisition has been obtained from. • Request ID. Identification of the client requesting the covering of a region. If a set of acquisitions have been generated for a client that wants to cover a region, this field has the same ID as the client in all the acquisitions generated, although the acquisitions can cover other different zones. •Priority. Attribute associated to the importance of the region covered, that is, priority p. •Start time. Time elapsed since the considered epoch at which the acquisition begins. •Stop time. Time elapsed since the considered epoch at which the acquisition stops. • Status. Binary attribute, it takes value 1 if the acquisition is part of the planning solution, 0 otherwise. Lastly, the acquisitions considered have a minimum and maximum duration due to time-store capabilities of the satellite. In terms of minimum time, considering problems as the heating of the sensors or scheduling slew manoeuvrers, an acquisition should last, at least, 60 seconds and, as a maximum size, each of them must not exceed the 3600 seconds. Subregions Consider a region R and a set of acquisitions A={A1,A2,...,An} . These acquisitions do not necessarily cover all the area of R . The subregions set S={S1,S2,...,Ss} is all the divisions inside Rgenerated by A, that is, S={S⊆R:Sis a subregion obtained by Ai∈A}, ∀S∈S,S⊆R, ∃S∈S:S⊆A. (2.3) As an example of the generation of subregions, consider Figure 2.1. The region R (blue rectangle) is divided into multiple subregions using the acquisitions in red (Figure 2.1a). As previously said, the acquisition set does not need to cover all the region. The subregions generated are inside R , as shown in Figure 2.1b. Following the definition in Equation 2.3, all the subregions generated are not inside the set of the acquisitions. 2.2 Formulation 7 (a) (b) Figure 2.1 Region and acquisitions (a) and subregions generated by the acquisitions inside the region (b). 2.2 Formulation Let R={R1,R2,...,Rr} be a set of uncovered regions with different priorities, pi . Let A= {A1,A2,...,An} be a set of different acquisitions generated by a set of satellites C={C1,C2,...,Cc} with modes M={M1,M2,...Mm} defined with a roll ( φ ) and a pitch ( θ ) angle, Mi={(φ,θ)|φ,θ∈ [−π 4,π 4]} ∀Mi∈M . A has been obtained within a time interval T= [Ti,Tf] . The goal is to cover the maximum area possible in R taking into account the priorities pi , so that Ri has preference in the scheduling plan over Rjif pi>pj∀i,j=1,...,r MILP formulation In order to build a Mixed Integer Linear Programming (MILP) formulation for our problem, the subregion set S={S1,S2,...,Ss} is obtained with A and R following the definition of subregion in Equation 2.3. The parameter σi defines the area of Si and the parameter ρi defines the different priorities of each subregion, so that, ρi=pjif Si⊆Rj∀Si∈S∀Rj∈R.(2.4) In order to solve the problem, two parameters are introduced: Qi j =(1,if Siis covered by Aj 0,otherwise. ∀Si∈S,∀Aj∈A(2.5) Qmode k j =(1,if Akis incompatible with Aj 0,otherwise. ∀(Ak,Aj)∈A(2.6) Furthermore, two decision variables are defined: xk=(1,if Akis part of the solution 0,otherwise. ∀Ak∈A(2.7) yi=(1,if Siis covered by the solution 0,otherwise. ∀Si∈S(2.8) 8Chapter 2. Problem Setting The formulation of the problem presented below is tailored to the needs of Deimos, with the objective of covering the largest possible area of the regions with the highest priority, regardless of the region’s overall size. To achieve this, an iterative procedure is introduced to find the optimal solution. 1. The subregions are sorted by priority, from the highest to the lowest value of ρi. 2. Subregions with the same priority are selected and stored in the set S∗ , where the number of subregions in this set is s∗≤s. 3. For this selected set of subregions, the following formulation is used: max s∗ ∑ i=1 σiyi(2.9) s.t. Qmode k j xj≤1si xk=1,k=1,...,n,j=1,...,n.(2.10) yi=1if Qi jxj≥1,i=1,...,s∗,j=1,...,n.(2.11) The objective function (2.9) maximizes the sum of the areas of the subregions with the same priority that are covered by the selected acquisitions. Constraint (2.10) ensures that each acquisition in the solution set is compatible with the others. Constraint (2.11) ensures that a subregion is included in the solution if one or more acquisitions cover that subregion. 4. The acquisitions selected as part of the solution are fixed in subsequent iterations by adding new constraints to (2.10) and (2.11) . This means that in later iterations, the variable xk will represent the acquisitions already chosen. 5. Steps 3 and 4 are repeated with the next set of subregions of lower priority, continuing until all subregions are covered, or no more acquisitions can be added to the solution. In the constraints presented in (2.10) and (2.11) , the variable is part of the restriction domain, which cannot be introduced in standard MILP solvers. This is why we reformulate the problem (2.9)-(2.11) using big-M constraints:. max s∗ ∑ i=1 σiyi(2.12) s.t. Qmode k j xj≤M(1−xk),k=1,...,n,j=1,...,n.(2.13) Myi≥Qi jxj−1+ε,i=1,...,s∗,j=1,...,n.(2.14) M(1−yi)≥1−Qi jxj,i=1,...,s∗,j=1,...,n.(2.15) Constraint (2.13) replaces constraint (2.10) and constraints (2.14) and (2.15) replace constraint (2.11) . The constant M is a large number, typically chosen to be significantly larger than any possible value of the decision variables, and it is related to a dimension of the problem. It needs to be large enough to satisfy the original constraints, but it should also be as small as possible to avoid numerical instability issues in the solver. For example, M could be chosen to be the sum of the sizes of sets S and A . Constant ε represents a tolerance related to the precision of the integers used in the solver. It ensures that the solver can handle the constraints properly while still achieving an optimal solution. Typically, ε is chosen to be a small positive value within the range (0,1) . This constant is important to ensure yi=1in the cases where Qi jxj=1in (2.14) and (2.15). Subregions calculation and matrices Parameter Q allows describing the problem in a MILP format to be later implemented into a solver. This matrix is obtained simultaneously with the calculation of the subregions. This calculation is 2.2 Formulation 9 performed using the computational geometry Elementary Polygons, integrated in MatLab. Let Sk be a subregion of S and Aj an acquisition of A . According to (2.3) , the two boolean operations that define a new set of subregions using Ajare: 1. Intersection. Let I={I1,I2,...,Ii} be the set of subregions generated by the intersection Sk∩Aj . This means that each subregion I in I satisfies the condition I∩Aj=I , indicating that each I∈Iis entirely contained within Aj. 2. Subtraction. Let B={B1,B2,...,Bb} be the set of subregions generated by the operation Sk−I . This implies that each subregion B in B satisfies the condition B∩Aj=/0 , meaning that each B∈Blies entirely outside of Aj. Figure 2.2 shows the boolean operations described in the previous points. If this process is (a) (b) Figure 2.2 Intersection and subtraction of a subregion and an acquisition (a) generating Iand B. repeated ∀Sk∈S (from now on, Sold ), a new set of subregions (from now on, Snew ) and a new matrix Q are obtained. The calculation of a new set of subregions along with the calculation of Q is presented in Algorithm 1. To obtain all the subregions that A generates in R , Algorithm 2 is presented. This algorithm computes the subregions generated in R and for every acquisition in A . As we will discuss later, this part of the process is a computational bottleneck, since the number of computations grows exponentially. 10 Chapter 2. Problem Setting Algorithm 1 Subregion generation with a single acquisition Require: Sold,Qold,Aj Snew ←− /0 Qnew ←− /0 for Sk∈Sold do I=Sk∩Aj B=Sk−I for I∈Ido Snew ←−I for each A∈Ado ▷Copy attributes for all A Qnew(I,A) = Qold(Sk,A) end for Qnew(I,Aj) = 1 end for for B∈Bdo Snew ←−B for each A∈Ado ▷Copy attributes for all A Qnew(B,A) = Qold(Sk,A) end for end for end for Algorithm 2 Subregions calculation Require: R,A Q(size(R),size(A)) = 0 S←−R for a∈Ado (S,Q) = SRandQ(S,Q,a)▷This is Algorithm 1 end for Lastly, the matrix Qmode k j must be obtained. Let Ak and Aj represent two different acquisitions from A that start at ts and end at tf . Let φk and φj be the roll angles, and θk and θj the pitch angles of the satellite at the moment the acquisitions are obtained. Let ˙ φ and ˙ θ be the angular velocities of the satellite for changing the roll and pitch directions, and let tpoint be the time required to adjust the camera when the satellite changes its attitude or captures another acquisition. The condition that ensures the compatibility between Akand Aj, that is, Qmode k j =0, is given by max{tf k,tf j}−min{ts k,ts j}≥(tf k−ts k)+(tf j−ts j)+tslew k,j,(2.16) where tslew k,jis calculated as tslew k,j=max |φk−φj| ˙ φ,|θk−θj| ˙ θ+tpoint.(2.17) 3 Input Data In this chapter, we introduce the Earth Observation Satellite Scheduling Problem (EOSSP) input data. The problem involves scheduling a set of acquisitions while adhering to temporal constraints, making the process of obtaining these acquisitions one of the key aspects. Besides acquisitions, the problem also depends on the ground track of the sensors onboard the satellites, also known as the swath. Therefore, information related to the orbital elements, the characteristics of onboard sensors, and the attitude of the satellites is necessary for resolving the EOSSP. All these elements, essential for solving the EOSSP, are depicted in Figure 3.1. This figure illustrates both the data acquisition required to solve the problem (blocks 1 and 2) and the problem resolution itself (block 3). It is important to highlight this division because Deimos Space already owns tools related to orbit propagation, attitude calculation, and sensor information (block 1), as well as tools related to computational geometry for obtaining acquisitions (block 2). However, Deimos lacks an algorithm to schedule these acquisitions. Therefore, the objective of this thesis is to develop a solver (block 3) capable of scheduling the acquisitions in the best possible way. Figure 3.1 Flowchart of the program used to solve the EOSSP. 3.1 Swath Generator The purpose of this block is to determine the ground track and swath by taking into account the sensor characteristics and the satellite’s attitude within a specified time interval. As previously mentioned, this block represents the fields of orbital mechanics and attitude dynamics within the scope of the work. The inputs and outputs of this block are detailed as follows. Inputs •Initial Epoch: The initial date on which the satellite starts to acquire data, denoted as ti. •Final Epoch: The final date on which the satellite finishes acquiring data, denoted as tf. •Time Step (τ): The minimum level of discretization for the problem. Combining the initial and final epochs with the time step, the satellite’s ground track can be discretized into N steps: N=⌊tf−ti τ+1⌋(3.1) 11 18 Chapter 4. Algorithms Algorithm 3 Greedy Algorithm Require: R(set of regions), A(set of acquisitions) Solution ←/0 Sort Rby priority for r∈Rdo r←r−∪(Solution) while A=/0 and r=/0 do A′←possibleAcq(A,Solution) if A′=/0 then break while end if Ar←r∩A′ if Ar=/0 then break while end if a←biggestArea(Ar) Solution ←Solution∪{a} A←A−{a} r←r−(r∩a) end while end for 4.2 GRASP Algorithm 19 4.2 GRASP Algorithm GRASP, or Greedy Randomized Adaptive Search Procedure, is an iterative multi-start metaheuristic algorithm used for finding good solutions to combinatorial optimization problems. Each iteration consists of two phases: a constructive phase, where a feasible solution is built in a greedy randomized manner, and a local search phase, where this solution is iteratively improved by exploring its neighbourhood. The best solution found over all iterations is then considered the final solution. As seen in Algorithm 4, the GRASP implemented for our problem has three different phases: constructive, repairing mechanism and local search (note that we need a third phase, the repairing mechanism, to adjust the solution obtained in the constructive phase). Algorithm 4 GRASP Algorithm Solution ←− Apply constructive phase to R(set of regions) and A(set of acquisitions) currentSolution ←−Solution Solution ←− Apply repairing mechanism to currentSolution currentSolution ←−Solution Solution ←− Apply local search to currentSolution Constructive Phase In this phase, an initial solution is constructed similarly to the Greedy Algorithm. An acquisition selected as part of the solution is chosen randomly from a Restricted Candidate List (RCL), which depends on the parameter α∈[0,1] . The closer α is to one, the larger the RCL. The construction of the RCL is denoted as RCL =CL(α) , where CL denotes the Candidate List, which is the full set of all possible acquisitions available for selection. The RCL, or Restricted Candidate List, is then a filtered subset of the Candidate List that includes only the best candidates according to a predefined criterion based on the parameter α. The process begins by sorting acquisitions based on the area they cover. The list generated (the Candidate List, CL ) is then reduced by filtering out certain candidates to form the RCL. The higher the value of α , the more acquisitions will be included in the RCL. From this resulting list, one acquisition is randomly selected as the new component of the solution and removed from the set of available acquisitions. As in the Greedy Algorithm, this process is repeated until there are no more acquisitions to be selected, no more acquisitions can be chosen due to operational constraints, or there are no more regions to cover. The key aspect of this process is the repetition of the described sequence. When no more acquisitions can be selected for the reasons mentioned, we have a complete solution, which is compared to previous solutions obtained using this process. The number of solutions, Nsol , represents the number of times the process is repeated, and the best solution found, according to a comparison criterion, is selected as the final solution for this phase. The parameter α has two extreme values. When α=1 , the RCL becomes equal to the CL, meaning the selection process is entirely random. In this case, sorting the candidate list by area is unnecessary since an acquisition will be chosen randomly from all available options. On the other hand, when α=0, we have the Greedy Algorithm. For the comparison method, all the covered areas of the regions are compared. If an improvement is shown in the region with the highest priority, the new solution is selected; otherwise, it is discarded. If the covered area does not change, the next region in priority queue is compared. This process of comparing old areas and new ones is denoted as comparisonMethod(areasold,areasnew) , and it is 20 Chapter 4. Algorithms detailed in Algorithm 5. This method of comparing areas by priority, and not the total area, must be done due to the objectives of Deimos Space presented in Chapter 2. Algorithm 5 Comparison Method Require: areasold,areasnew (both with the same length and sorted by priority) for i=1to length(areasold)do if areasnew(i)>areasold(i)then return true else if areasnew(i) = areasold(i)then continue else return false end if end for Algorithm 6 Constructive Phase Require: A,R,Nsol Sort Rby priority bestAreas ←0 while Nsol >0do Nsol ←Nsol −1 Solution ←/0 for r∈Rdo r←r−Solution A′←possibleAcq(A,Solution) while r=/0 and A′=/0 do Ar←r∩A′ if Ar=/0 then break while end if CL ←sort Arby area RCL ←CL(α) a←random candidate from RCL Solution ←Solution ∪{a} A′←A′−{a} r←r−(r∩a) end while end for areas ←calculateAreas(Solution) if comparisonMethod(bestAreas, areas) then bestAreas ←areas bestSolution ←Solution end if end while 4.2 GRASP Algorithm 21 Repairing mechanism The solution obtained during the constructive phase may not be optimal in terms of the number of acquisitions used. This inefficiency arises because some acquisitions might redundantly cover the same parts of regions due to the iterative updating of the regions. To address this issue, a repairing mechanism is introduced, as described in Algorithm 7, specifically an Erase Mechanism repairing process. This repairing method removes acquisitions that redundantly cover the same area while taking into account the priorities of the regions. It is crucial to consider these priorities because some acquisitions may cover regions with higher priority while also covering regions with lower priority, and vice-versa. Algorithm 7 Repairing Mechanism Require: Solution eraseAcq ←/0 for s∈Solution do auxSolution ←Solution - eraseAcq - {s} newArea ←area(auxSolution) if newArea = area(Solution) then eraseAcq ←eraseAcq ∪{s} end if end for Solution ←Solution - eraseAcq Local Search The local search phase, as outlined in Algorithm 8, is designed to identify a better solution starting from an initial solution, S . The primary objective of this phase is to explore the neighbourhood of solutions around Sto improve upon it. In this phase, the first step is to construct a neighbourhood of solutions, denoted as V(S) . This neighbourhood consists of all potential solutions that can be obtained by making small changes to the initial solution S . These changes could involve swapping acquisitions, adjusting schedules, or making other minor modifications that are feasible within the constraints of the problem. Once the neighbourhood V(S) is established, the next step is to evaluate the solutions within this set. The goal is to identify an improved solution from this neighbourhood that better meets the criteria of the problem, such as reducing conflicts or enhancing coverage. The local search process then selects the best solution from V(S) and potentially updates S with this new solution if it offers a significant improvement. This process is iterative and continues until no further improvements can be found within the neighbourhood or a stopping criterion is met. This approach can avoid getting stuck in local maxima or minima, but it also increases the computational time. Algorithm 8 Local search Require: S (initial solution) Snew ←−Apply swapping to S Sbest ←−Apply adding to Snew Overall, the local search phase is focused on refining the initial solution by exploring neighbouring solutions and selecting the most promising ones to improve the quality of the acquisition schedules. Due to the randomized nature of the GRASP algorithm, this phase is crucial for discovering better 22 Chapter 4. Algorithms solutions. The previous theses introduced in the first chapter, such as [ 3 ] and [ 6 ], do not explore the repairing mechanism and the local search phases of GRASP algorithms. In this work, the local search phase is divided into two subproblems: 1. Swapping, as described in Algorithm 9, is a process where the set S represents the scheduled acquisitions, and N denotes the set of non-scheduled acquisitions. The neighbourhood V(S) is constructed by swapping an acquisition s∈S with a non-scheduled acquisition n∈N , provided that the priority of n is greater than or equal to the priority of s ( p(n)≥p(s) ), and n is compatible with the set S′=S−{s} . This phase of the local search involves replacing one of the scheduled acquisitions with one from the list of non-scheduled acquisitions to generate a new solution. It is important to note that non-scheduled acquisitions are explored for every scheduled acquisition removed from the solution. Consequently, the neighbourhood generated could be immense. In the worst-case scenario, if all the non-scheduled acquisitions could be included in the solution every time an acquisition from the solution set is swapped, the neighbourhood would have |S|×|N|possible solutions to compare. 2. Adding, as described in Algorithm 10, involves the set S of scheduled acquisitions and N of non-scheduled acquisitions. The neighbourhood V(S) is constructed by adding acquisitions from N to S . This local search method aims to increase the space of acquisitions scheduled. The primary difference compared to the swapping method is that now, for every acquisition included in the solution, the non-scheduled acquisitions set is explored to add more acquisitions to the set apart from the one included. When no more acquisitions can be added, the new set is included in the neighbourhood, and the process is repeated with the new acquisition from the list of non-scheduled acquisitions. In this case, the worst-case scenario is to include all the non-scheduled acquisitions in the set of scheduled acquisitions, resulting in |N| possible sets of acquisitions to compare in the neighbourhood. The increase in computation time is due to the necessity of comparing |N|non-scheduled acquisitions |N|times. Algorithm 9 Swapping Method Require: S (acquisitions in the solution), N (acquisitions outside the solution) V(S)←/0 for s∈Sdo S′←S−{s} for n∈Ndo if n is compatible with S’ and p(n) ≥p(s) then V(S)←V(S)∪{S′∪{n}} end if end for end for bestAreas ←areas(S) for v∈V(S)do if comparisonMethod(bestAreas, areas(v)) then bestAreas ←areas(v) Snew ←v end if end for 4.2 GRASP Algorithm 23 Algorithm 10 Adding Method Require: S (acquisitions in the solution), N (acquisitions outside the solution) V(S)←/0 N′←N for n∈Ndo N′←N′−{n} if n is compatible with S then S′←/0 S′←S∪{n} for n′∈N′do if n’ is compatible with S’ then S′←S′∪{n′} end if end for V(S)←V(S)∪S′ end if end for bestAreas ←areas(S) for v∈V(S)do if comparisonMethod(bestAreas, areas(v)) then bestAreas ←areas(v) Snew ←v end if end for 24 Chapter 4. Algorithms 4.3 Branch and Bound algorithm Branch and Bound is an iterative algorithmic strategy used for solving hard combinatorial optimization problems [ 4 ]. It involves dividing the problem into smaller subproblems, or branches, and systematically exploring the solution space to identify the optimal solution. The technique operates by recursively partitioning the search space into smaller regions, called branches, and maintaining bounds on the possible solutions within each branch. It utilizes a combination of depth-first or breadth-first search and pruning to efficiently explore the solution space while avoiding the need to examine all possible combinations. At each step, the algorithm evaluates partial solutions and uses lower and upper bounds to determine whether a branch can be pruned, thus avoiding the need for further exploration. This pruning process significantly reduces the computational effort required to find the optimal solution. The computation time order of the branch and bound algorithm can be estimated with three key components: 1. Branching Factor (b). This factor represents the average number of leafs states generated from each branch state. This depends on the problem and the branching strategy used. 2. Depth of the Solution (d), that represents the depth of the search tree from the root node to the deepest leaf node that contains a solution. 3. Bounding Function. This function determines whether a subtree can be pruned based on the current state and the best solution found so far. Thus, the time complexity is O(bd) . For our specific problem, each acquisition can or cannot be included in the solution, so b=2 and, in the worst-case scenario, the depth of the solution is equal to the number of acquisitions, that is, d=n . Then, the time complexity of the algorithm (excluding the bounding function) is O(2n) . To reduce the computation time, two pruning techniques are implemented in the algorithm: 1. Upper bound. Let S be the best solution found, with covered area A(S) . Let B be a branch generated and L the leafs the branch generates. The best scenario is that each leaf in L could be included in the solution, that is, S∗=B∪L . If A(S∗)≤A(S) , the branch B is pruned, as it cannot upgrade the best solution found. 2. Pruning factor ( β ). Let B be a branch and L be the leafs the branch generates. L is sorted in descending order of area covered, A(B∪Li)Li∈L . The set of leafs L′=L(β) are the new leafs selected as the leafs of B , with β∈(0,1] . The closer β to one, the closer L′ is to L . This factor is analogous to the α parameter in GRASP Algorithm 4, but now all acquisitions in the leaf can be considered as part of the new solutions, rather than selecting only one randomly. The bounding function is possibleAcq, explained in Algorithm 11, with a similar output as Algorithms 3 and 4 and, obviously, following the equation 2.16. For a branch B , the function generates the set of leafs L from the set of acquisitions A compatible with B , and then, returns the set of leafs L′=L(β) . This function also generates the upper bound explained previously. The branch and bound process is presented in Algorithms 12 and 13. To summarize the key characteristics of the algorithms, Table 4.1 provides a comparative overview highlighting their theoretical advantages and disadvantages. 4.3 Branch and Bound algorithm 25 Algorithm 11 Bounding function (possibleAcq) Require: branch, A(set of acquisitions), β(pruning factor) for a∈Ado if {a}compatible with branch then leaf ←leaf ∪{a} end if end for Sort each component in leaf by area covered leaf ←leaf(β) UpperBound ←area covered by (branch ∪leaf) Algorithm 12 Branch and Bound search (BnBSearch) Require: branch, leaf, A, UpperBound if leaf is empty then solution ←branch if area covered by branch is better than current area covered else if UpperBound ≤bestArea then leaf ←/0 end if for i∈leaf do newBranch ←branch∪{i} [newLeaf, newUpperBound] =possibleAcq(newBranch,A,β) BnBSearch(newBranch, newLeaf, A, newUpperBound) end for end if Algorithm 13 Branch and Bound process Require: R(set of regions) ,A(set of acquisitions) sort Rby priority solution ←/0 for r∈Rdo branch ←solution [leaf,UpperBound] = possibleAcq(branch,A,β) solution ←BnBSearch(branch, leaf, A, UpperBound) end for 26 Chapter 4. Algorithms Table 4.1 Advantages and disadvantages of the algorithms. Advantages Disadvantages Fast execution time Does not guarantee optimal solution Simple to implement Can get stuck in local optimaGreedy Works well for problems where local choices lead to good global solutions Suboptimal in complex or large problem spaces Combines greedy heuristics with randomness for better exploration More computationally demanding than Greedy due to multiple iterations and local search Flexible and adaptable to different problem instances Dependent on parameters for performanceGRASP Tends to produce better solutions than pure Greed No guarantee of finding the global optimal solution Guarantees an optimal solution in small problems Computationally intensive, especially for large problem spaces Systematically explores all possibilities Exponential time complexity in the worst-case scenarioBranch and Bound Efficient pruning reduces the solution space when possible Memory usage can grow significantly with the size of the problem 5 Experiments In this chapter, a set of experiments are conducted to test the capabilities of the different algorithms developed. The two main parameters to compare results is the computation time, Tcomp , and the Gap between the area covered by the algorithm, A(alg), and the area covered by the optimal solution, A(opt). The Gap definition in (5.1) is used to compare the total area covered and also the area of each region with different priority. Gap =A(opt)−A(alg) A(opt)∗100.(5.1) Additionally, since the type of satellites and sensors affects the solution, their specific characteristics are introduced at the beginning of the experiments. There are also other factors, such as the duration of acquisitions, that remain consistent across all experiments. For this thesis, and based on data provided by Deimos Space, it is established that the minimum acquisition duration is 60 seconds, while the maximum duration is 3600 seconds. The data from the different regions of the Earth to cover with the satellites have been downloaded from several database on internet such as Global Administrative Areas (GADM) [ 2 ] or Natural Earth [1]. 5.1 Experiment 1: testing the capabilities of the algorithms In this experiment, a set of clients wants to cover different areas in the Mediterranean Sea within a day. The objective is to cover the maximum area possible taking into account the priority levels between the regions. The areas to be covered with their priority are presented in Table 5.1. The higher the priority, the more importance of a region in the scheduling plan and, as can be seen, the priority is totally independent from the area of the regions. Table 5.1 Experiment 1: regions to be covered. Region Area [km2]Priority Greece 131957 1 Crete 8336 4 Cyprus 9250 3 Corsica 8680 5 Sardinia 2490 2 Balearic Islands 4992 6 27 34 Chapter 5. Experiments Computational time for = 0.5 10 50 100 200 500 Number of solutions (#) 0 200 400 600 800 1000 1200 1400 1600 Computational time (s) Constructive phase Repairing mechanism Local search Figure 5.11 Experiment 1: computation time for α=0.5and different number of iterations. 5.1 Experiment 1: testing the capabilities of the algorithms 35 Computational time for = 0.75 10 50 100 200 500 Number of solutions (#) 0 200 400 600 800 1000 1200 1400 Computational time (s) Constructive phase Repairing mechanism Local search Figure 5.12 Experiment 1: computation time for α=0.75 and different number of iterations. Branch and Bound The branch and bound results and the computational time are presented in Figures 5.13 and 5.14. The solutions obtained are better as the pruning factor increases, something that was expected. The interesting results in this type of algorithms is the computation time, shown in Figure 5.14. As β increases, the computation time increases in an order of magnitude. When comparing the increase in computation time with the improvement in the solution, the Gap does not decrease by orders of magnitude as β increases. The worst case is shown with the values β=0.01 and β=0.03 , in which only Greece and Crete has lower Gap, but there are not even the highest priority regions. 0 2 4 6 8 10 12 14 Gap [%] Balearic Islands Corsica Crete Cyprus Sardinia Greece Regions = 0.01 = 0.03 = 0.05 Figure 5.13 Experiment 1: Branch and Bound Gap for different β. 36 Chapter 5. Experiments 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Prunin factor ( ) 102 103 104 105 Computational time (s) Computational time fo Branch and Bound Algorithm Figure 5.14 Experiment 1: Computation time for Branch and Bound algorithm depending on the pruning factor. Comparison among different solutions To provide an overview of the results from this first experiment, the best set of solutions for each type of algorithm is compared in this section. The results are presented for the Greedy algorithm, the best value of Nsol for each α considered in the GRASP algorithm (denoted as GRASP(α,Nsol) ), and for the Branch and Bound algorithm with β=0.05 . The Gap values for these algorithms, along with their computation times, are shown in Table 5.3. Table 5.3 Experiment 1: comparison among the solutions obtained from different types of algorithms. Optimal Solution Greedy GRASP(0.25,50) BnB(0.05) Area Covered [%] Tcomp[s]Gap [%] Tcomp[s]Gap [%] Tcomp[s]Gap [%] Tcomp[s] Balearic Islands 95.28 5.89 0.55 5.89 Corsica 99.79 0.254 1.2 6.02 Crete 78.02 13.85 16.149 0.54 Cyprus 54.07 1.13 7.489 14.1 Sardinia 19.75 -0.131 2.08 2.08 Greece 46.24 4084.74 9.86 18.49 33.56 403.518 2.49 128867 Using the metric covering the largest possible area of the region with the highest priority, it is clear that the best solution is achieved using the GRASP method with α=0.25 and Nsol =50 , with a Gap of 0.55%. The problem with this solution is that the coverage of regions like Crete and Greece has high Gap values. Regarding computation time, it is 403 seconds, nearly 7 minutes, which is a good value. If the goal is to seek a balance of coverage across all areas, meaning that all customers are equally satisfied, the Branch and Bound solution with β=0.05 should be selected, where the Gap values are very similar to each other and reasonably low. The problem with this solution is the computation time, on the order of 105 seconds, which is unfeasible if the solution needs to be delivered to customers on the same day as the request. For a quick solution with good results, the best option is to choose the Greedy solution, which has acceptable Gap values and the lowest computation time of all, 18.49 seconds. 5.2 Experiment 2: Europe and USA coverage mission 37 5.2 Experiment 2: Europe and USA coverage mission The capabilities of the algorithms have been tested in the previous experiment, but it is not a real case scenario for Deimos Space or, at least, it has not been one yet. The satellite DEIMOS-1 has a swath width of 600km, so it is typically required for covering large areas, which is the purpose of the mission introduced in this section. In this experiment, DEIMOS-1 has the objective of covering the area of Europe and USA within 30 days as the primary mission and, as a secondary mission in the same timespan, it has to acquire other small regions such as islands or cities. The orbital elements in the epoch are the same as in Experiment 1 for SAT-1, but now the Angle of Vision of the sensor is 35ºand the roll and pitch angles are 0ºthroughout the entire mission. Table 5.4 Experiment 2: list of regions to be covered. Region Area [km2] Priority Europe 10,53 ∗1062 USA 9,8∗1062 Mediterranean Islands 1,6∗1051 Caribean Islands 2,3∗1051 Ceuta and Melilla 30 1 Madagascar 5,8∗1051 Figure 5.15 Experiment 2: regions to be covered. 38 Chapter 5. Experiments Deimos Space obtains the orbital elements of the satellite using a propagator within three days since the epoch, so, if the mission has a time horizon of 30 days, the algorithms will work 10 times to find the solution, once every 3 days. The number of acquisitions obtained are presented in Figure 5.16. As described in Chapter 3, the global program used to obtain all the elements to get the solution is lineal, so the acquisitions are obtained every time for the whole regions, taking into account also the parts that can be covered by previous acquisitions scheduled. 0 5 10 15 20 25 30 Day 80 81 82 83 84 85 86 87 88 Acquisitions Number of acquisitions obtained Figure 5.16 Experiment 2: number of acquisitions obtained every 3 days. 5.2 Experiment 2: Europe and USA coverage mission 39 Optimal solution The optimal solution for this problem is obtained 10 times as well as the solution of the algorithms, once every time a set of possible acquisitions are calculated. The subregions are obtained taking into account the parts of the regions that are already covered in previous days. The number of subregions obtained through the mission are shown in Figure 5.17. The number of subregions decays exponentially as more area of the region is covered. 3 6 9 12 15 18 21 24 27 30 Day 0 50 100 150 200 250 300 350 Subregions Number of subregions generated Figure 5.17 Experiment 2: number of subregions generated every 3 days. The area covered by the satellite is presented in Figure 5.18. Although the priority of Europe and USA is larger, the proportion of the area covered is lower than the area of the regions with priority 1, a solution expected due to the smaller area to be covered. The computation time is shown in Figure 5.19. The time to get the solution decays exponentially and then it stabilizes around the day 9. Again, the computation time share of the optimal solution is dominated by the algorithm obtaining the subregions and matrices. 40 Chapter 5. Experiments 3 6 9 12 15 18 21 24 27 30 Day 80 82 84 86 88 90 92 94 96 98 100 Area covered [%] Area covered by optimal solution P = 2 P = 1 Figure 5.18 Experiment 2: optimal coverage. Computation time optimal solution. Total time: 171.6497 s 5 10 15 20 25 30 Day 0 5 10 15 20 25 30 35 Time [s] Subregions Solver Figure 5.19 Experiment 2: computation time of the optimal solution. 5.2 Experiment 2: Europe and USA coverage mission 41 Greedy algorithm The Relative Percentage Deviation of the Greedy algorithm solution respect to the optimal solution is shown in Figure 5.20. As can be seen, the variation is smaller than 0.01%. The computation time of the algorithm, as shown in Figure 5.21, is lower than the computation time of the MILP in the first 9 days. Then the computational time of each day is larger than the optimal solution, with a peak in the day 24. 3 6 9 12 15 18 21 24 27 30 Day 0 1 2 3 4 5 6 7 8 9 10 Gap [%] 10-3 Gap Greedy algorithm P = 1 P = 2 Figure 5.20 Experiment 2: Greedy Gap. 3 6 9 12 15 18 21 24 27 30 Day 16 18 20 22 24 26 28 Time [s] Computation time Greedy algorithm. Total time: 202.4916 s Figure 5.21 Experiment 2: computation time of Greedy algorithm. 42 Chapter 5. Experiments GRASP algorithm The solutions for this experiment using the GRASP algorithm were obtained with values of α= 0.25,0.5, and 0.75 , and with the number of iterations in the constructive phase set to 10, 50, 100, and 200. A solution derived from a specific α and Nsol combination is denoted as GRASP(α,Nsol) and α . Figures 5.22, 5.23, and 5.24 display the Gap values for each day, corresponding to different values of Nsol . Similar to the previous experiment, the best results were obtained for α=0.25 , with the GRASP(0.25,100) solution having the lowest Gap value for the region where P=1 . A common observation across all graphs is that the Gap for the region with P=1 consistently tends towards zero, whereas the Gap for P=2 exhibits more fluctuations, with no clear trend. For α=0.5 and 0.75 , negative Gap values are observed during the initial days for P=2 , indicating that the GRASP algorithm found a better solution for the secondary mission, but at the cost of a worse solution for the primary mission. In the case of GRASP(0.5,200) , the secondary mission improves by 4% while the primary mission deteriorates by 1%. However, there are instances where the improvement in the P=2 mission is less significant than the degradation of the P=1 mission, as seen with the GRASP(0.75,50) solution. Overall, it appears that increasing the value of Nsol results in better outcomes for the P=2mission for a given value of α. 3 6 9 12 15 18 21 24 27 30 Day -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Gap [%] Gap GRASP(0.25,10) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Gap [%] Gap GRASP(0.25,50) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Gap [%] Gap GRASP(0.25,100) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day -0.5 0 0.5 1 1.5 2 2.5 3 3.5 Gap [%] Gap GRASP(0.25,200) P = 1 P = 2 Figure 5.22 Experiment 2: Gap for α=0.25 and Nsol =10,50,100 and 200. 5.2 Experiment 2: Europe and USA coverage mission 43 3 6 9 12 15 18 21 24 27 30 Day 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Gap [%] Gap GRASP(0.5,10) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day -1 -0.5 0 0.5 1 1.5 2 2.5 Gap [%] Gap GRASP(0.5,50) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day 0 0.5 1 1.5 2 2.5 3 3.5 4 Gap [%] Gap GRASP(0.5,100) P = 1 P = 2 3 6 9 12 15 18 21 24 27 30 Day -4 -3 -2 -1 0 1 2 Gap [%] Gap GRASP(0.5,200) P = 1 P = 2 Figure 5.23 Experiment 2: Gap for α=0.5and Nsol =10,50,100 and 200. 50 Chapter 5. Experiments 5.3 Experiment 3: coverage of small areas using a constellation The scope of this experiment is to test the capabilities of a constellation of satellites covering smaller areas with a backup mission of a bigger area. The goal is covering as much area as possible as the regions with different priorities are requested. The initial time horizon is three days from the epoch and then, increase the timespan day by day. If this is done 10 times, the total number of days of the mission is 12. As in the previous experiments, the acquisitions are obtained previously to save global computational time. The constellation is formed by 4 satellites with the orbit characteristics of DEIMOS-2 (described in Table 5.2 in SAT-2) and an Angle of Vision of 2º, that leads to a width of the swath of, closely, 12 km. The mean anomaly (M) of the satellites at the epoch is 0º,90º, 180ºand 270º. The acquisitions obtained in this experiment are descending, that is, are obtained from North to South. The pitch of the satellites is 0ºand the roll angle is between -45ºand 45º, with an angle step of 1º. The regions covered, their priority and the day they are requested are in Table 5.5. The priorities are selected according to the day of the request, having higher priority those regions requested first. The backup mission, or secondary mission, is the coverage of Africa, with priority 1. The rest of the areas are included in the primary mission, although they do not have the same priority. Table 5.5 Experiment 3: regions requested. Days accumulated New regions requested Priority Type of mission 3 Africa 1 Backup 3 Mediterranean Island, Ceuta and Melilla 100 Primary 4 Canary Islands 90 Primary 5 Caribbean Islands 80 Primary 6 Korean peninsula and Japan archipelago 70 Primary 7 Tasmania 60 Primary 8 Madagascar 50 Primary 9 Iceland 40 Primary 10 Baja California Peninsula 30 Primary 11 New Zeeland 20 Primary 12 Philippines 10 Primary Figure 5.32 shows the number of acquisitions obtained over the days. The initial three days witnessed the highest number of acquisitions, followed by a drop on the fourth day. Subsequently, the number of acquisitions increased as new regions were introduced. This pattern is expected as the purpose of the study is to test the algorithm’s capabilities across the entire region, including parts that could potentially be covered if they were incorporated into the algorithm. Optimal solution The optimal solution for this simulation has not been obtained due to the computation time employed to acquire the first iteration. The computer run the simulation for 4 days and the first iteration had not been completed, so there will be no comparison with the rest of the algorithms. Thus, the Gap parameter is not used as a comparison method. This leads to the vast difference in computation time of the optimal solution depending on the number of acquisitions considered. 5.3 Experiment 3: coverage of small areas using a constellation 51 3 4 5 6 7 8 9 10 11 12 Day 1000 1500 2000 2500 3000 3500 4000 4500 Number of Acquisitions Number of acquisitions obtained each day Figure 5.32 Experiment 3: number of acquisitions. Greedy algorithm The Figure 5.33 shows the area covered by the Greedy algorithm for the primary and backup missions. First of all, both primary and backup mission are not fully complete within 12 days. On the one hand, backup mission coverage follows a linear trend within the time simulated. On the other hand, the primary mission has more valleys and peaks, with a great drop in the coverage at day 5. This drop can be explained due to the inclusion of more area to be covered with lower priority whereas other areas with higher priority were still not fully covered. 3 4 5 6 7 8 9 10 11 12 Day 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered by Greedy algorithm Primary Backup Figure 5.33 Experiment 3: area covered by greedy algorithm for the primary and backup missions. The Figure 5.34 precisely shows the previously described behaviour in the primary mission. As progressively larger areas with lower priority are included, the remaining uncovered areas from previous regions become smaller, causing the growth to be slower compared to other experiments. 52 Chapter 5. Experiments It can be observed that the regions with higher priority are the ones that increase their covered area the fastest, with the exception of Tasmania (P = 60) and Iceland (P = 40), regions with significantly lower area than other regions considered. 3 4 5 6 7 8 9 10 11 12 Day 0 10 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered with different priorities by Greedy algorithm P = 100 P = 90 P = 80 P = 70 P = 60 P = 50 P = 40 P = 30 P = 20 P = 10 Figure 5.34 Experiment 3: area covered by the greedy algorithm for the primary mission showing the difference in priorities. Lastly, the computation time of the greedy algorithm is presented in Figure 5.35, along with the time limit of the simulation computation. If the algorithm is finding the solution for a propagated orbit over a day, the results cannot take more than a day to be obtained. As can be seen, the computation times each day rises as the number of acquisitions increases. 3 4 5 6 7 8 9 10 11 12 Day 103 104 105 Time [s] Computation time Greedy algorithm. Total time: 105161.3611 s Computation time Time limit Figure 5.35 Experiment 3: computation time of the greedy algorithm. 5.3 Experiment 3: coverage of small areas using a constellation 53 GRASP algorithm For the GRASP algorithm something similar to the optimal solution happened. The iterations needed to construct the iterative phase were running or 4 days and they were not completed, so there will be no solution for the GRASP algorithm to this experiment. Branch and Bound algorithm For the Branch and Bound algorithm, two different values of the pruning factor are selected to test the capabilities of the algorithm in this type of mission, one being one order of magnitude bigger than the other: β=0.0001 and β=0.001 . Despite the order followed showing the previous results, for this part the computation time is presented in Figure 5.36. It can be observed that for 3456789 Day 104 105 Time [s] Computation time BnB algorithm = 0.001 = 0.0001 Time limit Figure 5.36 Experiment 3: computation time of the branch and bound algorithm. both beta values, on day number 9, the simulation time limit is exceeded, meaning the algorithm would take more than a day to provide the solution intended for that specific day. On the other hand, it is interesting to note that for the first four days, for both beta values, the Branch and Bound algorithm returns a solution in less time than the one obtained with the greedy algorithm, before the computation time increases drastically. This is why the results only include information up to day nine, as it did not make sense to continue the simulation once the time limit had already been exceeded. The results are very similar to those obtained with the Greedy algorithm as can be observed in Figures 5.37 and 5.38, with the main drawback being that the complete simulation could not be finished. The figures practically show the same results but the solution with higher pruning factor takes more time to be obtained. 54 Chapter 5. Experiments 3456789 Day 10 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered by BnB algorithm, = 0.0001 Primary Backup 3456789 Day 10 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered by BnB algorithm, = 0.001 Primary Backup Figure 5.37 Experiment 3: area covered by branch and bound algorithm for the primary and backup missions. 3456789 Day 0 10 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered with different priorities, BnB with = 0.0001 P = 100 P = 90 P = 80 P = 70 P = 60 P = 50 P = 40 3456789 Day 0 10 20 30 40 50 60 70 80 90 100 Area covered [%] Area covered with different priorities, BnB with = 0.001 P = 100 P = 90 P = 80 P = 70 P = 60 P = 50 P = 40 Figure 5.38 Experiment 3: area covered by branch and bound algorithm for the primary mission showing the difference in priorities. 6 Conclusions and Future Work The main objective of this work was to develop a set of algorithms capable of scheduling acquisitions captured by a constellation of satellites, considering different pointing modes and constraints while maintaining solutions obtained that cannot be erase from the set of acquisitions scheduled. This objective has been achieved, as three approaches for solving the arising combinatorial optimization problem have been presented and tested in different types of missions. We observed in the experiments that none of the algorithms proposed is superior to the others for solving the EOSSP, as the algorithm performance strongly depends on the mission conditions, such as the number of satellites and the duration of the mission. The optimal solution can be found for scheduling acquisitions if the problem size, in terms of the number of acquisitions, is small enough. However, as the problem size increases, the computation time required to compute the optimal solution becomes impractical. The Greedy algorithm provides the best performance in terms of computation time, making it ideal for quickly generating a solution. The GRASP algorithm improves upon the solution obtained by the Greedy approach, offering better coverage in terms of area, though at the cost of increased computation time. Furthermore, GRASP’s two control parameters allow for greater flexibility in adjusting the quality of the solutions. The most important parameter in this algorithm is α , as smaller values of α tend to result in better solutions. However, as α increases, larger values of Nsol are required to achieve good results and, just like any other metaheuristic, no optimality guarantee is provided. Finally, the Branch and Bound algorithm has the highest computation time but can still be useful in specific situations. The pruning factor is a key parameter that makes this algorithm practical in real-world scenarios. One of the most important findings regarding this solution technique is the significant increase in computation time as the pruning factor β increases, with only marginal improvements in the solution quality for higher values of β. In conclusion, there is no clear winner among the algorithms presented. The specific characteristics of each mission will determine which of the techniques and algorithms is the most suitable for the task at hand. In terms of future work, several avenues could be explored to further improve the performance and flexibility of the algorithms developed. One possibility is to investigate the integration of machine learning techniques to predict and adapt the most suitable algorithm or combination of algorithms for specific mission characteristics, such as region size, priority, or satellite constellation parameters. Additionally, more advanced metaheuristic algorithms, like genetic algorithms or ant colony optimization, could be studied to enhance the efficiency of the scheduling process, especially in larger-scale missions. Another interesting direction would be the development of hybrid algorithms that combine the strengths of multiple methods (e.g., combining the speed of Greedy algorithms with the accuracy of Branch and Bound). Furthermore, it would be valuable to 55 56 Chapter 6. Conclusions and Future Work explore real-time scheduling capabilities, where the system can adapt to dynamic changes in satellite availability or mission objectives. Lastly, collaboration with satellite operators to implement and validate these algorithms in real-world scenarios would provide critical insights into their practical applicability and pave the way for further refinements based on operational data. List of Figures 2.1 Region and acquisitions (a) and subregions generated by the acquisitions inside the region (b) 7 2.2 Intersection and subtraction of a subregion and an acquisition (a) generating Iand B9 3.1 Flowchart of the program used to solve the EOSSP 11 3.2 Swath calculation flowchart 13 3.3 Example of acquisition generation: Spain and its influence region 14 3.4 Example of acquisition generation: preliminary acquisitions obtained to cover Spain 15 3.5 Example of acquisition generation: final acquisitions obtained to cover Spain 15 5.1 Experiment 1: acquisitions generated 28 5.2 Experiment 1: subregions generated by the acquisitions 29 5.3 Experiment 1: acquisitions scheduled by optimal solution 29 5.4 Experiment 1: area covered of each region by optimal solution 30 5.5 Experiment 1: acquisitions scheduled by greedy solution 30 5.6 Experiment 1: Gap of greedy solution 31 5.7 Experiment 1: GRASP Gap with α=0.25 and different number of iterations 32 5.8 Experiment 1: GRASP Gap with α=0.5and different number of iterations 32 5.9 Experiment 1: GRASP Gap with α=0.75 and different number of iterations 33 5.10 Experiment 1: computation time for α=0.25 and different number of iterations 33 5.11 Experiment 1: computation time for α=0.5and different number of iterations 34 5.12 Experiment 1: computation time for α=0.75 and different number of iterations 35 5.13 Experiment 1: Branch and Bound Gap for different β35 5.14 Experiment 1: Computation time for Branch and Bound algorithm depending on the pruning factor 36 5.15 Experiment 2: regions to be covered 37 5.16 Experiment 2: number of acquisitions obtained every 3 days 38 5.17 Experiment 2: number of subregions generated every 3 days 39 5.18 Experiment 2: optimal coverage 40 5.19 Experiment 2: computation time of the optimal solution 40 5.20 Experiment 2: Greedy Gap 41 5.21 Experiment 2: computation time of Greedy algorithm 41 5.22 Experiment 2: Gap for α=0.25 and Nsol =10,50,100 and 200 42 5.23 Experiment 2: Gap for α=0.5and Nsol =10,50,100 and 200 43 5.24 Experiment 2: Gap for α=0.75 and Nsol =10,50,100 and 200 44 5.25 Experiment 2: computation time for α=0.25 and Nsol =10,50,100 and 200 45 5.26 Experiment 2: computation time for α=0.5and Nsol =10,50,100 and 200 46 5.27 Experiment 2: computation time for α=0.75 and Nsol =10,50,100 and 200 47 57 58 List of Figures 5.28 Experiment 2: Branch and Bound Gap with β=0.01 47 5.29 Experiment 2: Branch and Bound Gap with β=0.05 48 5.30 Experiment 2: computation time Branch and Bound algorithm with β=0.01 48 5.31 Experiment 2: computation time Branch and Bound algorithm with β=0.05 49 5.32 Experiment 3: number of acquisitions 51 5.33 Experiment 3: area covered by greedy algorithm for the primary and backup missions 51 5.34 Experiment 3: area covered by the greedy algorithm for the primary mission showing the difference in priorities 52 5.35 Experiment 3: computation time of the greedy algorithm 52 5.36 Experiment 3: computation time of the branch and bound algorithm 53 5.37 Experiment 3: area covered by branch and bound algorithm for the primary and backup missions 54 5.38 Experiment 3: area covered by branch and bound algorithm for the primary mission showing the difference in priorities 54 List of Tables 4.1 Advantages and disadvantages of the algorithms 26 5.1 Experiment 1: regions to be covered 27 5.2 Orbital elements of SAT-1 and SAT-2 at epoch 22/5/24 00:00:00 UTM 28 5.3 Experiment 1: comparison among the solutions obtained from different types of algorithms 36 5.4 Experiment 2: list of regions to be covered 37 5.5 Experiment 3: regions requested 50 59