Quantum Optimization Methods for Satellite Mission Planning

Satellite mission planning for Earth observation satellites is a combinatorial optimization problem that consists of selecting the optimal subset of imaging requests, subject to constraints, to be fulfilled during an orbit pass of a satellite. The ever-growing amount of satellites in orbit underscores the need to operate them efficiently, which requires solving many instances of the problem in short periods of time. However, current classical algorithms often fail to find the global optimum or take too long to execute. Here, we approach the problem from a quantum computing point of view, which offers a promising alternative that could lead to significant improvements in solution quality or execution speed in the future. To this end, we study a planning problem with a variety of intricate constraints and discuss methods to encode them for quantum computers. Additionally, we experimentally assess the performance of quantum annealing and the quantum approximate optimization algorithm on a realistic and diverse dataset. Our results identify key aspects like graph connectivity and constraint structure that influence the performance of the methods. We explore the limits of today’s quantum algorithms and hardware, providing bounds on the problems that can be currently solved successfully and showing how the solution degrades as the complexity grows. This work aims to serve as a baseline for further research in the field and establish realistic expectations on current quantum optimization capabilities.


I. INTRODUCTION
T HE Satellite Mission Planning Problem (SMPP) is a crit- ical issue in the aerospace sector [1].Satellite operators face the task of determining the optimal subset of images to be captured during a satellite's orbital pass, based on a set of client requests and a defined value metric.The challenge is compounded by constraints such as geographical proximity incompatibilities, onboard storage limitations, and specific image configuration requirements, which make it unfeasible to collect all images.Fig. 1 depicts a simple scenario of a SMPP where the satellite has to choose from the set of all requested images (rectangles) which ones to take (green) and which to leave unattended (red).The problem's complexity increases when considering the continuous influx of requests, thereby necessitating frequent plan updates and repeated execution of the planning algorithm within short time frames.This requirement underscores the importance of execution speed.Current industry challenges such as extending the problem to satellite constellations, incorporating additional constraints like cloud coverage, or real-time onboard re-planning [1], increase computational demands even further.
The SMPP and its extensions exhibit high combinatorial complexity, as even its simpler versions can be cast to variants of the knapsack problem [2], which is known to be NPcomplete [3].This poses significant challenges for classical computing methods, even for moderately sized instances with few constraints.Traditional classical approaches have employed exact integer linear programming algorithms based on branch-and-bound/cut methods [2], [4], [5], (meta-)heuristic algorithms [6]- [9], and machine learning based techniques [10].The latter are typically favored due to their ability to handle large problems and comply with execution-time constraints.
Quantum computing (QC) is emerging as a promising alternative for addressing hard optimization problems by leveraging properties of quantum physics [11], [12].Several approaches are currently being investigated by the community, of which quantum annealing (QA) [13] and variational quantum algorithms such as the quantum approximate optimization algorithm (QAOA) [14] are among the most prominent.These methods, coupled with the appropriate hardware, may provide near-to mid-term benefits such as substantial computational speedups, better-quality solutions or a reduction in energy consumption.
In this paper, we study the SMPP from the QC point of view, using approaches from two of the leading and available paradigms: quantum annealers and gate-based devices.To do that, we present a formulation in which a variety of hard constraints are analyzed, paying extra attention to the more challenging ones.We use 31 different problem instances to evaluate our formulation as well as the performance of quantum annealing and QAOA methods on currently available quantum hardware and simulators.Six of these instances were directly sourced from the well-known SPOT5 dataset [2], while the remaining 25 were generated based on SPOT5 using an instance reductor specifically implemented for this study.
The rest of the paper is structured as follows.Section II provides a comprehensive background related to the problem under investigation.The mathematical models developed to efficiently formulate the problem from the classical and quantum points of view are presented in Section III.Section IV describes the quantum algorithms evaluated and the criteria for performance assessment and comparison.Section V details the dataset, experimental design, conducted tests, and a discussion of the findings.The paper culminates with Section VI, which draws conclusions and outlines potential avenues for future research.

II. BACKGROUND
For years, quantum technologies have had a significant influence in the space sector [15], with many resources dedicated to secure satellite communications and QKD [16], [17] or quantum sensing [18].Notably, there has been a recent broadening of focus within the industry towards quantum machine learning and optimization [19].This has led several major space agencies to establish initiatives like the ESA's Quantum Computing for Earth Observation1 (QC4EO), NASA's Quantum Artificial Intelligence Laboratory2 (QuAIL) or DLR's Quantum Computing Initiative3 (QCI) to propel advancements in research and development within this domain.
Due to the importance of the SMPP in the aerospace industry and the limitations in terms of computational speed and solution quality when classical computers are used to solve it, methods from quantum computing appear as an attractive alternative.However, most of the published work does not deal with the SMPP per se, but with other related topics such as satellite routing, debris removal, or constellation optimization [20], [21].We take the scarcity of literature as motivation for this research and devote the rest of this section to analyzing the published related work.
The seminal work in this field was published in 2020 by Stollenwerk et al. [22], [23], which deals for the first time with the SMPP using the quantum annealing paradigm.They investigate the potential and maturity of then-current quantum computers to solve real-world problems by carrying out an experimental study on a reduced number of small instances of the SMPP.
Another notable contribution is [24], where the authors demonstrate potential for quantum advantage by using a hybridized quantum-enhanced reinforcement learning algorithm and compare it to greedy and classical optimization algorithms on a satellite mission planning problem.On another note, in [25] variational quantum algorithms such as the variational quantum eigensolver [26] and QAOA are used to solve several simple instances of the SMPP with and without simulated noise.
Additionally, the authors of the present paper published a preliminary study in [27], which analyzed the problem formulation efficiency together with the performance of the state-of-the-art quantum annealing solvers at the time.
Thus, this work builds on that preliminary research and improves it by studying a larger set of problems with more diverse constraints.The study includes a theoretical complexity analysis of the models and a more in-depth experimental study using gate-based and annealing methods.More concretely, the main elements of this paper with respect to the state of the art and our preliminary findings can be summarized as follows: • The dataset employed in our research, taken from realistic simulations of a satellite mission, has been extended to cover a wide range of problem sizes and constraint types for a more in-depth study.• We perform a detailed result analysis, measuring the quality of the solutions via the approximation ratio against reference classical solutions on the actual optimization function value, rather than using proxies such as the energy or only the best solution.Lastly, it is noteworthy to acknowledge the existence of publications such as [28], [29], which employ quantum-inspired genetic algorithms.These methodologies belong to the field of quantum-inspired evolutionary computation, which are classical algorithms augmented by principles derived from quantum physics for their design and thus are outside the scope of this paper.

III. MODELLING
In this section, we present the problem and the mathematical formulations used in this paper.First, in Section III-A we introduce the key aspects of the problem we are going to solve.Then, in Section III-B, we present a classical mathematical model for it.Subsequently, for each type of constraint we encounter in our model, we discuss in Section III-C possible penalty encodings and indicate their efficiency in terms of additional variables and quadratic terms.Finally, using these derivations we arrive at our formulation for the experiments which will be conducted in Section V.

A. PROBLEM DESCRIPTION
In this paper, we work with data from the SPOT5 satellite, which has three cameras on board and can satisfy two types of requests: mono images, for which only one of the cameras must be used, and stereo images, for which two of the cameras are required to produce the image.Each request has a weight that represents its value, and it may also have a capacity weight if it needs to be stored on the satellite's disk.The capacity value can vary depending on the camera used to take the request.This definition is subject to the following constraints: a) Constraints enforcing that each request can only be collected once.b) Binary constraints that represent incompatibilities due to geographical proximity among certain pairs (p, q), where p and q are pairs of (request, camera).We define, for each instance of the dataset, a set S 2 containing these forbidden pairs.c) Ternary constraints that represent the instantaneous data flow restrictions of the instruments, prohibiting taking more than 2 out of 3 requests at once for certain triples (p, q, r), where p, q and r are pairs of (request, camera).We define, for each instance of the dataset, a set S 3 containing these forbidden triples.d) N-ary (capacity) constraints that represent the fact that the satellite has a limited amount of disk space available to store images before relaying them back to earth.This type of constraint is present only in some of the instances.For a more in-depth description of the dataset, the reader can refer to [2].

B. CLASSICAL MODEL
The classical formulation of the SMPP used in this paper builds upon the ideas provided on our previous research [27], and it is stated using mathematical programming as follows: Let x i,j be the binary decision variables, defined as: where i ∈ {1, 2, . . ., N } is the index representing the image requests, N being the total amount of requests and j ∈ {1, 2, 3, 4} the identifier of the camera.There are three physical cameras that can take mono images and we define camera 4 to represent the combined use of cameras 1 and 3 to take stereo images.The necessary variables are created based on the specific data instance, following the notation mentioned above.The objective function to be optimized can then be defined as: where w i ∈ R ≥0 is the weight or value of fulfilling request i.Note that although our task is to maximize the value, we can express it as the minimization of the negative of the objective function.This optimization is subject to the following constraints: where constraints (3a)-(3d) refer to the four constraints mentioned above, respectively.Here C ∈ Z ≥0 is the total capacity of the disk and c i,j ∈ R ≥0 the amount of space a request i uses when taken with camera j.For some (i, j) we effectively have c i,j = 0 due to the corresponding orbit points being close enough to a ground station for immediate relay, hence no storage requirement.Importantly, this makes constraint (3d) take quite different forms for different instances.As a final remark, note that we can always flatten our decision variables into a one dimensional vector through the map x i,j → x i ′ , which is what we will use in practice for our implementation of the experiments.

C. EFFICIENT QUANTUM MODEL
Quantum optimization methods usually require to express the problem in a quadratic unconstrained binary optimization (QUBO) formulation due to their connection to the Ising model [30] and adiabatic quantum computing [31], which serves as inspiration for many current methods like QA and QAOA.However, mapping a classical optimization problem to a QUBO is, in general, a difficult task that the quantum optimization community is commonly faced with [32].Additionally, the efficiency of this mapping is critical to the performance of current quantum computers, given their limitations in scale, connectivity, and stability.To this end, we will focus our attention on the encoding of constraints as penalties.In a nutshell, penalties are a method of encoding constraints in the objective function by adding a large M ∈ R multiplied by an at most quadratic polynomial that represents the constraint.Thus, on a minimization problem, adding a large positive M makes the objective function larger when the constraint is violated.
We analyze the efficiency of the encodings by looking at the two factors that are within our control and impact the problem size the most: the number of slack variables (ancillary qubits) and the number of quadratic interaction terms (qubit couplings) introduced by the encoding methods.
Constraints (3a) can be trivially converted to penalties by taking the sum of the two-by-two products of the variables involved in each constraint, resulting in no variable overhead and k where k is the number of summation terms.In our case, k ≤ 3 since the image is either stereo (only possible to take the image with camera 4) or mono (possible to take the image with at most three different cameras).Therefore, at most three interaction terms are introduced for each constraint.We proceed in the same fashion with constraints (3b), as they have a similar structure.In this case, k = 2 and thus no additional terms are introduced and one interaction term is added for each constraint.Constraint (3c) can be efficiently treated in a custom manner.First, we write an equivalent cubic penalty term M x p x q x r and then reduce it to quadratic following Boros et al. [33,Sec 4.4].A single slack variable s 1 replaces the pair x q x r , and one extra interaction term is added, resulting in the following quadratic penalty, M x p s 1 + M (x q x r − 2x q s 1 − 2x r s 1 + 3s 1 ). (4) This method has advantages for constraints such as (3c) because it involves fewer terms than the more general binary expansion (discussed next) and introduces a single slack variable per constraint.Furthermore, if the pair x q x r is present in multiple constraints, it can be replaced by the same slack variable in all of them, preventing the introduction of additional slack variables.
To deal with constraint (3d), which, as mentioned above, can take different forms depending on the data, we resort to the standard binary expansion.For ease of reading, we take p = (i, j) and rewrite the constraint as: where we sum over every pair p for which c p ̸ = 0, with P such pairs in total.Now we transform this inequality into an equality.To this end, we add binary slack variables s d such that the binary representation reaches or exceeds C, that is, we find the minimum D such that This results in an overhead of D = ⌈log 2 C⌉ slack variables.The constraint can now be written as: which allows the use of the general penalty method: Therefore, with this method we are adding O(log C) ancilla qubits and O (P + log C) 2 quadratic interactions for each constraint of that type.As a note, although constraint (3c) is a special case of (5), the well-defined form of the former (three variables and fixed coefficients) compared to the more variable nature of the latter compels us to treat them in quite different ways.While the Boros approach could in principle also be used to address constraint (3d), two significant challenges arise.First, we would need to construct the equivalent polynomial to be reduced, which in general is a highly combinatorial problem.And even so, faced with a simple case in which c p = 1 ∀p, where the polynomial to be penalized is simply the sum of all products composed of C +1 distinct x p variables, the number of new slack variables added scales with O(P C ).This makes the method unsuitable for widespread use in generic constraints.

IV. METHODS
In this section, we briefly describe the algorithms and solvers we use, as well as the measure of solution quality we chose for this study.Sections IV-A and IV-B outline the quantum annealing process and the QAOA algorithm, respectively.In Section IV-C we describe the characteristics of D-Wave's hybrid solver, and finally, in Section IV-D we describe the criteria we use to evaluate solution quality.
A. QUANTUM ANNEALING Quantum annealing processors are especially suitable for combinatorial optimization problems.They are thought to exploit quantum tunneling, entanglement, and superposition to escape local minima and explore a large solution space [34].In QA, a quantum system is made to evolve under a Hamiltonian that interpolates between a simple Hamiltonian called mixer Hamiltonian, H M , of which the ground state is known, and an Ising Hamiltonian that encodes the desired problem, called cost Hamiltonian H C .A usual choice for the mixer Hamiltonian is H M = i X i , where X i is the x Pauli operator on the i-th qubit, and |+ . . .+⟩ is its ground state.On the other hand, the Ising Hamiltonian has the following form: where Z i is the z Pauli operator on the i-th qubit, h i is its on-site energy and J i,j the coupling between qubits i, j.Therefore, we have the following total Hamiltonian: where f (t) ∈ [0, 1] is the interpolation function, usually defined as a polynomial function of time t.By means of the adiabatic theorem [35], if we start at the ground state of H M and evolve the system slowly over time to H C , it will remain at the ground state and thus provide a solution to the minimization problem.
In QA, the choice of interpolation function aims to correctly retrieve the ground state of H C in the end, but not necessarily traversing the ground state of H (f (t)) in intermediate steps, which allows for faster evolution.
We can minimize the corresponding function on classical binary variables x i (obtained by transforming Z i → 1 − 2x i ).This is by analogy called energy, and can be written as: where Q is a matrix that can be assumed symmetric without loss of generality.A QUBO consists in the minimization of such functions.Current commercial quantum annealers, such as the ones provided by D-Wave, natively solve this kind of problems.The limitation to polynomials of order two comes from the fact that the Hamiltonian in (8) only couples qubits two by two, and this is intrinsic to the hardware.Additionally, the hardware is not able to directly couple every pair of qubits, so an important feature of the quantum processing unit (QPU) is its topology, i.e. the configuration of qubit couplings.Due to this connectivity limitation, most problems require complex embedding algorithms to represent them on the QPU.Usually, several physical qubits are required to encode a logical one, creating chains that have to stay coherent, which is a challenge for current devices.For this study, two different QPUs have been used: Advantage_system6.4 (Advantage) and the recently released Advantage2_prototype2.2 (Advantage2).The topologies of these devices are Pegasus and Zephyr, respectively, the latter having higher connectivity [36], [37].

B. QUANTUM APPROXIMATE OPTIMIZATION ALGORITHM
The quantum approximate optimization algorithm [14] is a hybrid algorithm that combines gate-based quantum computing with classical variational-parameter optimization in order to solve combinatorial optimization problems.
This algorithm can be understood as a discretized approximation of adiabatic evolution, provided the right choice of parameters.It does so by constructing a variational ansatz that alternates ℓ times the unitary operators encoding the cost Hamiltonian H C and the mixer H M (same as in QA), which results in: where the 2ℓ variational parameters are γ = (γ 1 , γ 2 , . . ., γ ℓ ), and β = (β 1 , β 2 , . . ., β ℓ ), with γ j ∈ [0, 2π) and This anstaz is applied to an initial state |ψ 0 ⟩, which is the ground state of the mixer H M and is generally easy to obtain.Thus we have the state: We denote E ℓ as the expected value of the Hamiltonian H C when acting on this quantum state.This expectation is given by: This quantity is minimized via a classical optimization routine.That is, we aim to find the optimal parameters: The solution is then obtained by measuring the state |Ψ ℓ (γ * , β * )⟩ in the computational basis.In principle, this should allow to obtain an energy equal to or close to the ground state energy of H C .The positive integer ℓ controls how discrete the evolution is.As ℓ → ∞, the closer to adiabatic the evolution can be, and, in general, the better the results.However, as ℓ increases, so does the number of parameters and the circuit depth, making optimization more difficult and time consuming.In order to compensate for this problem, we use the parameter-fixing strategy introduced in [38].Starting from a ℓ = 1 version of QAOA with random initialization (which works for shallow circuits), the 2ℓ parameters optimized for a ℓ-layer QAOA are used as initialization values of a (ℓ + 1)-layer QAOA along with γ ℓ+1 = β ℓ+1 = 0.

C. LEAP BQM HYBRID SOLVER
The Leap Binary Quadratic Model Hybrid (LeapBQMHybrid) is a component of the Hybrid Solver Service (HSS, [39]) offered by D-Wave, which is a collection of algorithms created by this company.In order to solve optimization problems that quantum processors are unable to directly handle, these hybrid approaches imbricate classical and quantum processing [40].As of this writing, the solvers included in HSS allow the addressing of Binary Quadratic Models (BQM), Discrete Quadratic Models (DQM) and Constrained Quadratic Models (CQM).The reason for using LeapBQMHybrid over the other available solvers is that the problem dealt with in this paper is defined entirely in binary variables.
Going deeper into the algorithm, the LeapBQMHybrid workflow is split into two distinct phases.To be more exact, the BQM formulation is first introduced as input into a classical front end.The solver then runs a predetermined number of parallel processing threads, each of which consists of a quantum module (QM) and a heuristic module (HM).On the one hand, HM is devoted to tackling the problem by means of traditional state-of-the-art heuristic techniques.On the other hand, QM helps HM by executing a variety of quantum queries with the objective of guiding the search towards promising regions of the search space.Eventually, QM can also improve the existing solutions through these queries.
Finally, LeapBQMHybrid provides the user with the optimal solution found among the group of generated threads.
It is interesting to mention that LeapBQMHybrid resorts to the latest D-Wave quantum computer.At the time this research was conducted, this device was the Advantage_system6.4,which is made up of 5616 qubits organized in a Pegasus topology.We depict in Fig. 2 the main workflow of this solver, which is based in the work published in [41].Readers interested in further details about the LeapBQMHybrid might peruse the D-Wave related report [40].

D. SOLUTION QUALITY ASSESSMENT CRITERIA
Benchmarking the performance of quantum optimization algorithms and hardware is still an open research question without a clear answer [42].In this study, we have adopted a similar approach to prior research in the field, using the approximation ratio (AR) as a metric for how good a solution vector x is.We define AR as follows: where a solution vector is deemed feasible if it respects the constraints in (3), and F max denotes the maximum objective function value achievable by any feasible solution.
It is important to note that from the runs of the algorithms, we get (n + s)-long vectors, where the first n components are the actual decision variables and the remaining s are slack terms introduced by the penalties.Therefore, we define the solutions x to be the first n components of these vectors.This definition allows us to evaluate the quality of solution vectors in the context of the optimization problem at hand, hence independently of the specifics of the QUBO encoding utilized, providing a standardized measure to assess their performance.

V. EXPERIMENTATION AND RESULTS
This section is devoted to the description of the main elements of the experiments conducted and the results obtained.In Section V-A we introduce the characteristics of the dataset utilized in our study and how we generated it.Then, in Section V-B we detail our design choices for the experiments performed.Finally, in Section V-C we discuss and analyze the results.Fig. 3 illustrates the implemented experimental pipeline.

A. DATA AND INSTANCE REDUCTOR
We have developed a Python script for the automatic production of SMPP instances as part of the current study.The implemented mechanism, coined SPOTReductor, works as follows: first, an existing instance of the problem (in SPOT5 format) is introduced to SPOTReductor together with the number of requests desired for the new instance.After that, SPOTReductor reduces the input problem by randomly selecting constraints and keeping the requests involved in them until the desired number of requests is reached.Then, the maximum capacity of the satellite is established as half of the capacity that would be used if all requests were taken.Finally, to generate the instances without capacity constraints, the capacity requirements for each request are simply set to zero and the maximum capacity is removed.This way, we can study more in-depth the influence of the capacities in the problem.
In all, the testing planned for this study has 31 distinct instances.Of them, 6 have been obtained from the SPOT5 benchmark dataset and an additional 25 instances have been created through SPOTReductor.With these instances, we cover a range of jobs from 3 to 209, however, it is worth noting that it is not trivial to order the instances by difficulty, as their intrinsic variety in terms of request type, amount of ternary constraints, capacities, etc. play a vital role in this regard.Here we chose to order them by the amount of requests for readability.Additionally, to improve this study's reproducibility, all generated cases along with the obtained results are openly available in [43].

B. EXPERIMENTAL SETTING
To conduct the experiments, we have used our QUBO/Ising formulation as described in Section III-C on the 31 instances detailed in Table 1.We obtained reference solutions (and hence F max ) from exactly solving the classical formulation with First, we generate the dataset by combining original SPOT5 instances with reduced ones.Then, we formulate the classical model, translate it to our proposed efficient QUBO encoding and obtain the equivalent Ising model.We get the reference solutions with OR-Tools and run the quantum/hybrid subroutines a total of 5 times per model and instance.Finally, we compute the approximation ratio for the best found solutions as well as the average from the 2000 samples of each run executed.

TABLE 1.
Main characteristics of the used instances, ordered by increasing number of requests.For each instance, we depict the number of total and stereo requests, the amount of total and ternary constraints, and the number of independent and interaction terms in our formulation.The naming convention is as follows: if an instance is from the original dataset [2], we keep its original name, which consists of a series of digits (e.g. 8, 404, 1502).If the instance was generated by our reductor, we use the following pattern: gDDD(c) (e.g.g003, g003c), where the character g is used to flag that the instance was generated, then three digits represent the amount of requests in the instance, and finally, a character c is present for instances that have a capacity constraint.Google OR-Tools [44].To build the QUBO matrices, we chose to set the value of all penalty coefficients in each instance to M = i w i + 1, which is a standard choice used in the literature if no prior information about the instances is known [32].To account for the probabilistic nature of the methods, we run each combination of instance and solver 5 times.Finally, we report the expected value of the AR of the solutions and the AR of the best sampled solution.To make QAOA comparable to QA in terms of the best solution, we sample from the QAOA statevector as many times as shots run for QA.

ID
For the experiments with QA, we tested three different D-Wave solvers: Advantage, Advantage2, and LeapBQMHybrid.For the two QPU-only based solvers, we have left all parameters at their default values except for the number of reads, which we have set to 2000.For the hybrid solver, only the time_limit parameter is tunable, which we chose to set to the default value computed by D-Wave for each instance.
For QAOA, starting with ℓ = 1 and using the parameterfixing strategy described in section IV-B, we run the algorithm for ℓ = 1, . . ., 10 for each instance.For the mixer Hamiltonian, we employ the standard H M = i X i , and thus the initial state is set to be its ground state, |ψ 0 ⟩ = |+ . . .+⟩.We use COBYLA as the optimizer [45], configured with a tolerance of 10 −6 as stopping criteria.We report results for ℓ = 1 and ℓ = 10.
To select the initial parameters from which to start the optimization process for ℓ = 1, we randomly pick for each instance 5 pairs of (γ 1 , β 1 ), run QAOA with depth ℓ = 1 for each of them and keep the best performing pair among the five, ranked by the approximation ratio of the output.
The executions with QAOA have been carried out in an ideal statevector simulator, namely Pennylane's default_qubit device [46].We refrained from using real hardware in the QAOA experiments as the number of requests that would need to be made to the quantum computers would be prohibitive in terms of wait times and costs.

C. RESULTS
Let us now analyze the outputs from our experiments, which are summarized in Fig. 4. First, we remark that our classical exact optimizer achieved the optimal solutions for all instances and therefore we are not going to focus on the comparison of the quantum methods against it, but rather use these solutions as a reference to compare the quantum (hybrid) methods among themselves.
Beginning with LeapBQMHybrid, we see that it performs remarkably well across the board, achieving perfect approximation ratios for all but 4 of the considered instances.However, as described in Section IV-C, it is not clear how much of the processing happens on the QPU and what is its contribution.For this reason, we focus our analysis on the other methods.
In the case of the purely quantum annealing results, it is interesting to see that the Advantage2, while being still in development, performs better or comparably to the established Advantage QPU for the smaller instances.This suggests that there is an evolution in the quality of the processors and sheds a positive outlook on the problems we will be able to tackle when it is released in full scale.
As seen in the lower half of the figure, the Advantage solver is able to find the optimal solution for 13 of the non-capacity instances and 7 of the capacity ones, while Advantage2 is slightly less performant, with 12 and 5, respectively.However, the general trend with Advantage2 is that it often samples better solutions overall, as seen in the upper part of the figure, providing better average approximation ratios.
It is also noticeable that finding an embedding of a problem on the QA's QPUs does not guarantee that a reasonable solution will be found.This is evident from the larger capacity instances, where we see that from instance g025c onwards, even if the problem is sufficiently small to be embedded in the hardware, the approximation ratio is 0. However, this behaviour is expected: if a problem's connectivity is close to the limit of the QA hardware, resulting embedding is very complex, necessitating very long chains of physical qubits to represent a logical variable.These long chains are more unstable and it is very challenging to keep them coherent, thus resulting in errors that ruin the solution.
Another interesting behaviour is observed in the expected AR of capacity instances with QA, where we see jumps in performance even for problems with similar amounts of requests.This is due to the fact that not all requests have nonzero capacity, which is made evident in Fig. 5, where we can see the contribution of the capacity constraint by checking the difference between the capacity and non capacity versions of the instances.We see that g004c is heavily connected, meaning most images need to be stored on the satellite, while g005c is sparse since many images can be directly relayed and thus do not add extra connections.These results indicate that QA is sensitive to the sparsity of the graph.While it might seem that QAOA is not affected by this phenomenon, we stress that we run ideal simulations, and higher connectivity still poses a problem for real gate-based quantum computers.Usually, they do not have all their qubits coupled among themselves, thus necessitating a large overhead of operations to conform to the connectivity of the problem, increasing the depth of the resulting circuits and consequently being subject to more noise induced errors.In general, the sparsity of the graph is a factor that needs to be taken into account, especially for noisy intermediate-scale quantum (NISQ) [47] computers and algorithms.
Turning our attention to QAOA, we see that the 10 layer version performs better than 1 layer across all instances when expected AR is considered, and on all but one instance (g006) when comparing the best-solution AR.The reason a shallower QAOA might perform better can be attributed to three factors.Firstly, the stochastic nature of the sampling process.Secondly, the fact that we optimize the parameters with respect to the expected value, not the best solution.And thirdly, the parameter optimization algorithm is naturally more prone to getting trapped at local minima in higher dimensional search spaces, resulting in suboptimal parameters for circuits with more layers.
Remarkably, looking at the bottom row of Fig 4, we see that QAOA performs significantly better on instances with capacities than without, which means that potentially, the graph structure influences the quality of the solutions given by the algorithm, being more suitable for complete or close to complete graphs.This finding supports the evidence from [48], where the authors show that QAOA performs better on complete and regular graphs than on random ones.This fact might become useful in the long-term, when (if) fault tolerant computers become available, as on NISQ devices the noise mostly negates this feature.
We further propose that this phenomenon may be attributed to the fact that capacity-constrained problems naturally impose tighter constraints, resulting in a reduced number of feasible solutions.During the optimization process executed by QAOA, the algorithm tends to redistribute probability mass from nonfeasible states, characterized by notably penalized energies, to feasible states.Given the scarcity of feasible states in capacityconstrained problems, there exists a higher probability, albeit fortuitously, that the redistributed probability mass aligns with the globally optimal state.

VI. CONCLUSIONS AND FURTHER WORK
In this paper, we have experimentally assessed the performance of quantum annealing and gate-based optimization algorithms on an efficient formulation of the SMPP with intricate constraints and structure.We have resorted to a realistic benchmark dataset and shown how the quality of the solutions degrades with problem size, imposing practical limits on the instances that can currently be solved effectively and providing a comprehensive overview of the challenges that arise when dealing with the SMPP in particular.This work assesses the readiness of quantum optimization methods and hardware for the SMPP, focusing on the gate-based and annealing paradigms.
We have shown that it is possible to achieve good results (best solution AR ≥ 0.9) for moderately sized instances of relatively complex SMPPs using quantum optimization methods.With QA, we were able to tackle instances with sizes of up to about 80 requests (≈ 200 variables) when no capacities are considered and up to 10 requests (≈ 30 J J J J J J J J J J J $SSUR[LPDWLRQUDWLR :LWKRXWFDSDFLW\ /HDS%40+\EULG $GYDQWDJH $GYDQWDJH 4$2$OD\HU 4$2$OD\HUV J F J F J F J F J F J F J F J F J F J F J F J F J F J F ,QVWDQFH :LWKFDSDFLW\ J J J J J J J J J J J ,QVWDQFH $SSUR[LPDWLRQUDWLREHVW J F J F J F J F J F J F J F J F J F J F J F J F J F J F ,QVWDQFH FIGURE 4. Results for all solvers and instances considered.On the left side of the figure, we show the results for the instances without capacity constraints, while on the right side, the instances with capacities are depicted.On the top row, we represent the expected value of AR, while on the bottom row, the AR of the best solution is shown.All plots have 95% confidence interval error bars for the 5 runs of each experiment.Note that we did not include the classical solver's results, as their approximation ratio is always 1 without variance. variables) when capacities are included.Meanwhile, with QAOA, the sizes of the instances we are able to deal with are more modest, achieving reasonable results for instances of up to 6-8 requests (≈ 20 variables), both with and without capacities, although we were limited by the number of qubits we were able to simulate.
On another note, the new generation of D-Wave quantum annealers is showing promising progress, which leaves an optimistic outlook on future developments.With respect to the QAOA, even in ideal simulated environments, the performance is generally not on par with QA.However, there is evidence [49], [50] that suggests that, compared to QA, QAOA can be more robust to vanishing spectral gaps, and it is possible that adding more layers to the circuit or finding better optimization strategies could improve the results.
Additionally, it is important to remark that the parametrization greatly influences the quality of both methods.This leads us to believe that an exhaustive tuning of the parameters of the problem (value of the penalties, variable encodings, etc.) and models (number of layers in QAOA, optimizers, annealing time in QA, etc.) might improve the results significantly, but this poses significant challenges, not only computational but also with hardware availability, especially if the tuning needs to be done on today's scarcely available quantum computers.
Lastly, future research could be focused on extending the planning to multiple satellites that can cooperate, which increases the complexity and presents a more realistic scenario of the industry's challenges.Additionally, it would also be interesting to compare our results with other variational algorithms or run the QAOA with simulated noise models or even on real hardware.Basic questions such as determining whether quantum phenomena positively affect optimization are also yet to be explored.Finally, a growing field of research is quantum-inspired optimization, where other algorithms such g004 g004c g005 g005c FIGURE 5. QUBO matrices of instances g004(c) and g005(c) represented as graphs.Some capacity constraints introduce more complexity to the problem than others depending on whether or not the requests take up capacity (due to being immediately relayed or not).
as tensor networks could be explored.

FIGURE 1 .
FIGURE 1. SMPP diagram.A satellite is orbiting the Earth and has to choose which image requests, represented by rectangles, to capture.One possible solution to the planning problem is to take the green ones and discard or delay for the next orbital pass the red ones.

FIGURE 3 .
FIGURE 3. Schematic diagram of the pipeline implemented to conduct the experiments.First, we generate the dataset by combining original SPOT5 instances with reduced ones.Then, we formulate the classical model, translate it to our proposed efficient QUBO encoding and obtain the equivalent Ising model.We get the reference solutions with OR-Tools and run the quantum/hybrid subroutines a total of 5 times per model and instance.Finally, we compute the approximation ratio for the best found solutions as well as the average from the 2000 samples of each run executed.
We analytically compare the complexity of encoding the constraints via several methods and propose a final, efficient formulation for the problem.•We perform an experimental study on the performance of several approaches, namely, two purely quantum QA solvers, a quantum-classical hybrid method based on QA, and a gate-based QAOA that involves a classical optimization stage.
• Several types of challenging constraints are addressed in the study.Notably n-ary (capacity) and ternary (•