Cutting-Edge Trajectory Optimization through Quantum Annealing

: This paper introduces an innovative approach to explore the capabilities of Quantum Annealing (QA) for trajectory optimization in dynamic systems. The proposed method involves transforming trajectory optimization problems into equivalent binary optimization problems using Quadratic Unconstrained Binary Optimization (QUBO) representation. The procedure is general and adaptable, making it applicable to a wide range of optimal control problems that entail the satisfaction of dynamic, boundary, and path constraints. Speciﬁcally, the trajectory is discretized and approximated using polynomials. In contrast to the conventional approach of determining the polynomial degree ( n ) solely based on the number of boundary conditions, a speciﬁc factor is introduced in our method to augment the polynomial degree. As a result, the ultimate polynomial degree is calculated as a composite of two components: n = l + ( m − 1), where m denotes the count of boundary conditions and l signiﬁes the number of independent variables. By leveraging inverse dynamics, the control required to follow the approximated trajectory can be determined as a linear function of independent variables l . As a result, the optimization function, which is represented by the integral of the square of the control, can be formulated as a QUBO problem and the QA is employed to ﬁnd the optimal binary solutions.


Introduction
Trajectory optimization plays a pivotal role across various domains, including space exploration and satellite missions, where finding the most efficient trajectories is essential for mission success, fuel conservation, and resource optimization.Recent advancements in optimization techniques have been propelled by emerging artificial intelligence technologies, notably Genetic Algorithms (GA) and Neural Networks (NN).
In the field of satellite constellations, Savitri et al. employed GA in [1] to optimize trajectories, maximizing coverage while minimizing revisit time.Rughani et al. in [2] utilized GA to optimize orbital trajectories for spacecraft swarms, facilitating collaborative tasks in space and collision avoidance.Particle Swarm Optimization demonstrated its prowess in [3], addressing real-time guidance for autonomous lunar landings and hazard avoidance, with promising experimental results.In [4], Particle Swarm Optimization was exploited to maximize asteroid surface coverage using hopping trajectories in lowgravity environments.Application to asteroids Itokawa and Bennu showed significant surface coverage improvement.A comprehensive overview of genetic algorithms was presented in [5], providing insights into established algorithms, their implementations, and associated pros and cons.Neural Network techniques have also made significant contributions.Federici et al. applied reinforcement learning in [6] to autonomously guide a spacecraft during a mission to impact a binary asteroid, utilizing neural networks for optical observation processing and achieving a high success rate despite uncertainties.Scorsoglio et al. in [7] leveraged reinforcement meta-learning, combining hazard detection and avoidance with image and radar data to ensure safe landing points.Addressing challenges in asteroid exploration, Jiang et al. in [8] proposed the use of hopping rovers and explored deep reinforcement learning for 3D terrain path planning.Gaudet et al. in [9] focused on optimizing an asteroid hovering controller through reinforcement metalearning, enabling stable positioning even with less precise asteroid data.Combination of these methods can yield advantageous outcomes, as seen in [10,11], where frameworks synergized particle swarm optimization for trajectory optimization and extreme learning machines for precise gravitational acceleration approximation from asteroids.
The advancements in this field have unlocked the potential for more advanced and efficient trajectory optimization methods, emphasizing the role of computational efforts in realizing these breakthroughs.As a result, the importance of computer processing power has become paramount in the field of trajectory optimization.Especially the computational limitations have become a significant challenge in handling complex trajectory design problems.Traditional computational resources struggle to handle the massive calculations required for optimal trajectory design, especially for scenarios involving intricate navigation through varying gravitational fields and orbital constraints.
Quantum Computing (QC) has demonstrated the ability to outperform classical and more commonly used optimizers in specific cases [12].Some problems where QC provides an advantage over classical algorithms are factorization of prime numbers [13], unstructured database search [14],linear systems of equations [15], and simulation of quantum systems [16][17][18].Research is active, with many algorithms proposed in many different areas, such as differential equations [19], artificial intelligence and machine learning [20,21] and optimization [22,23].
While the aforementioned literature focuses primarily on gate-based quantum computing, which operates by applying a series of gates to manipulate qubits, similar to classical computers, another notable quantum computing model is QA.QA exploits and guides the natural evolution of a quantum system to solve combinatorial optimization problems.During QA, the system gradually transitions from an initial state to a low-energy state, similar to the cooling process in metallurgy, hence the term "annealing".QA promises more efficient optimization solutions, especially for large-scale and complex problems [24,25].In particular, current quantum annealers can handle computations with about an order of magnitude higher number of qubits than gate-based quantum machines.Proposed applications of QA include factorization of prime numbers [26], systems of polynomial equations [27], and differential equations [28,29].
The research is also exploring the application of quality control to practical engineering challenges, despite the limitations of current quantum annealers.In aerospace engineering, some applications have been proposed, such as addressing some artificial intelligence problems framed as combinatorial tasks [30], optimizing aircraft trajectories to prevent conflicts [31] and allowing agile planning of the Earth observation maximization of image acquisition [32].
The goal of this study is to explore the cutting-edge capabilities of QA and its potential applications.Although QA has garnered considerable attention in the optimization community, its potential in the specific context of space trajectory optimization has yet to be fully explored.The paper presents a transcription procedure that converts trajectory optimization problems into equivalent binary optimization problems using the QUBO representation.This versatile procedure can be applied to a wide range of optimal control problems, encompassing dynamic, boundary, and path constraints.Specifically, the trajectory is discretized and approximated using polynomials for each component.By employing inverse dynamics, the required control to follow the approximated trajectory can be determined.Consequently, the optimization function becomes the system's energy, represented by the integral of the square of the control.As this function depends on the form of the polynomial, the focus lies on optimizing the polynomial itself.This novel approach demonstrates significant potential for efficiently exploring extensive solution spaces, making it particularly promising for addressing complex trajectory design challenges in fields like space exploration and satellite missions.This paper is organized as follows.With the purpose of providing the appropriate context and establishing the necessary nomenclature, Section 2 provides a brief description of the basic principle behind Quantum Annealig.Section 3 presents an overview of the transcription procedure for trajectory optimization, elucidating the mathematical formulation of the problem, including the cost function and various types of constraints.Moving on to Section 4, the lunar landing and rendezvous problems and their simplifications are reported.The results of the transcription procedure for this specific problem are discussed in Section 5, with an emphasis on specific considerations and limitations.Finally, in Section 6, conclusions and potential future research directions in the field of QA for trajectory optimization are presented.

Quantum Annealing Optimization
The fundamental principle behind QA is to encode the optimization problem into the energy states of a quantum system.For this purpose, let us briefly recall some fundamental concepts of quantum physics: a quantum system can have discrete levels of energy, among which the one with minimum energy is called the ground state.The energy of a quantum system is represented by its Hamiltonian operator H, which is closely related to the time evolution of the system through the famous Schrödinger Equation [33].It is also possible to write a so-called eigenvalues problem of the Hamiltonian, as shown in Equation (1); this equation is sometimes called the Time Independent Schrödinger Equation.|ψ n (written here in the Dirac notation |• , used to indicate quantum states) are called eigenvectors (or eigenstates) of H, and E n are the eigenvalues that correspond to the energy values that the system can assume for n = {0, 1, 2, . ..}.Therefore, an optimization problem encoded into the Hamiltonian of a quantum system has the possible values of the cost function corresponding to the energy values E n of the Hamiltonian, with E 0 , the ground state, being the global minimum of the cost function.
Quantum annealers are typically represented by a network of qubits, for which the Hamiltonian assumes the form of a matrix.Qubits are the quantum equivalent of classical bits, but unlike classical bits, which can only represent zero or one, qubits can exist in the superposition of both states simultaneously.This allows QA exploration of multiple possible solutions to the optimization problem simultaneously.The optimization process in QA involves the adiabatic theorem from quantum mechanics, which states that a system in the ground state, whose Hamiltonian is slowly time-varying, remains in the ground state of the instantaneous Hamiltonian.Therefore, for every instant of time t from zero to annealing time T, the following equation holds: The optimization procedure starts by preparing the system in the ground state of a simple Hamiltonian, called driver Hamiltonian (H D ); the optimization problem is instead encoded in the energy eigenvalues of a problem Hamiltonian, H P ; then, the system is evolved adiabatically from H D to H P [30] in such a way that the system state at the end of the evolution is the ground state of the problem Hamiltonian.At this point, to obtain the solution, we need to only operate a measurement of the qubtis [34].Equation (3) shows a generic transition between the Dirver and Problem Hamiltonian, where A(t/T) and B(t/T) are arbitrary functions such that A(0) = 1, B(0) = 0 and A(1) = 0, B(1) = 1.In Equation ( 4), the state of the system is shown by the eigenvalues equation at the start and at the end of the annealing process, where ψ 0,D stands for the ground state of the Driver Hamiltonian and ψ 0,P is the ground state of the Problem Hamiltonian, the solution of the optimization problem.
To be solved by QA, an optimization problem must have a precise structure given by the Ising or QUBO formulations [35], which are described below.

•
Ising Formulation revolves around determining the minimum energy configuration of a system of interacting spins.Let us consider an Ising model with n spins denoted as σ = (σ 1 , σ 2 , . . ., σ n ), where each spin can be in one of two states: +1 (up) or −1 (down).The energy of the Ising model is expressed by the following equation: In Equation ( 5), a i represents the linear coupling coefficients for individual spins, while b ij represents the coupling coefficients for pairs of spins.The first sum accounts for the energy associated with the individual spin states, whereas the second sum reflects the energy due to interactions between spins.• QUBO Formulation is a powerful representation of optimization problems using a quadratic function of binary variables.We consider a problem with n binary variables denoted as x = (x 1 , x 2 , . . ., x n ), where each x i can take either zero or one.The objective is to minimize the following quadratic function: In Equation ( 6), x T represents the transpose of the binary variable vector x, and Q is an n × n matrix of coefficients.The goal is to find the assignment of binary values to x that minimizes quadratic function J q (x).This formulation is widely used in QA to transform complex optimization problems into a format compatible with quantum systems.
The Ising and QUBO formulations are closely related through a simple transformation that maps Ising spins (σ i ) to binary variables (x i ).In the Ising Formulation, the linear coupling coefficients (a i ) are transformed into the diagonal components of matrix Q in the QUBO Formulation, while the coupling coefficients for pairs of spins (b ij ) are represented by the off-diagonal components.It is important to recognize that these formulations are mathematical abstractions designed to conveniently express optimization problems in a general format.When implementing an optimization problem on a real quantum annealer, the possible mismatch between the logical lattice of the problem (how the variables mathematically interact with each other) and the physical connectivity pattern (how the qubits are actually connected with each other in the annealer) must be taken into account.In fact, quantum annealers are built with a specific graph topology that dictates which qubits are physically connected with other qubits.If the logical graph of the problem does not match the physical graph of the machine, a procedure called embedding must be employed to adapt the first to the second.The embedding works by creating chains of qubits in the physical graph that are forced to behave as a single variable, always assuming the same value; through these chains, it is possible to connect qubits not physically connected in the annealer topology and crafting an embedded optimization problem equivalent to the starting one.
As an extremely simple schematic to illustrate the concept, in Figure 1,a logical graph is shown where three variables (represented as circles) interact as follows: Variable 1 interacts with Variables 3 and 4. On the other hand, the physical graph, on the right, does not provide a connection between Qubits 1 and 4. Therefore, a possible embedding is to create a chain including Qubits 2 and 4, such that Connection 1-4 can now be implemented exploiting Connections 1-2 and 2-4.We notice that, although the logical problem includes three variables, the embedded problem must use four qubits; therefore, the embedding procedure comes at the cost of increasing the number of qubits needed to tackle the problem.
As a consequence, a problem that requires embedding could require an annealer with a number of qubits much larger than the number of variables.Moreover, utilizing longer chains in the embedding can increase susceptibility to noise and errors.Longer chains may become fragile, leading to potential disruptions that adversely affect the quality and reliability of the solution.Addressing these challenges is critical to unlocking the full potential of QA for solving complex optimization problems and realizing the transformative power of quantum computing.
One company that has made significant contributions in this area is D-Wave Systems, which is a leading company specializing in QA.Their quantum processors, known as quantum processing units (QPUs), are designed to perform QA computations with fully quantum systems (for more details, please refers to the following link: https://www.dwavesys.com/media/htjclcey/advantage_datasheet_v10.pdf,accessed on 6 November 2023).D-Wave's QPUs are used in various research and commercial applications to tackle realworld optimization challenges.D-Wave solvers are constantly being improved to address the limitations of quantum computing, and future advancements will likely include more sophisticated algorithms and improved processing capabilities to enable more efficient optimization.

Transcription Procedure
The goal of trajectory optimization is to find the optimal control inputs (u(t)) that steer a system to the desired final state (x f ) while considering the trade-offs between achieving the objective and minimizing the associated cost.Thus, given an initial condition (x 0 ) and a control input trajectory defined over finite interval t ∈ [t 0 , t f ], the cost function (J) in the Bolza form can be computed as follows: where χ(x(t 0 ), x(t f )) is the cost component depending on initial state x(t 0 ) and final state x(t f ) of the system, and , u(t), t)dt represents the integral of the instantaneous cost function L over the entire trajectory, which depends on state x(t), control input u(t), and time t.Moreover, the trajectory optimization problem must align with the underlying dynamics of the system.In the proposed method, the dynamic systems considered are linear or linearized; thus, the trajectory optimization problem can be expressed in the following comprehensive form: Equation ( 8) provides a concise representation of the trajectory optimization problem, encapsulating the goal of reaching a desired final state while considering the cost incurred along the trajectory, subject to the system's dynamic behavior and initial condition.In particular, x = [s(t), ṡ(t)] T represents the state vector, resulting in a total of six components.This vector is composed of both position (s(t) = [x, y, z] T ) and velocity (ṡ(t) = [v x , v y , v z ] T ) along the three trajectory components.
The proposed transcription procedure aims to map the common trajectory optimization problem (Equation ( 8)) into a binary optimization problem in a QUBO form (Equation ( 6)), which can be effectively addressed using QA techniques.Central to the transcription procedure is the discretization of the trajectory by representing the state variable and its derivatives through polynomial functions across finite interval t ∈ [t 0 , t f ].
In contrast to the usual practice of determining the polynomial degree (n) based on the number of boundary conditions, our method increases the polynomial degree by a specific factor.Consequently, the ultimate polynomial degree is obtained as the combination of two components: n = l + (m − 1), where m represents the number of boundary conditions and l signifies the count of independent variables.These independent variables are essentially the coefficients of the polynomials that play a crucial role, acting as pivotal crossing points capable of altering the trajectory's shape.This core concept forms the foundation of the proposed optimization approach, where both the trajectory and its derivatives can be expressed as functions of l independent variables.In our procedure, we chose to minimize the number of variables, thus selecting l = 1, leading to polynomial degree n being simply equal to m.The decision to set l = 1 was primarily dictated by a desire to keep the approach simple and practical, making it easier to convert to the QUBO format.Moreover, it is worth noting that in other works, a similar choice was made to set "L" to 1 [3,36].This approach is commonly adopted in trajectory optimization research for its practical advantages.Thus, the main goal of the article is not to add complexity to the underlying mathematics, but rather to test QA with entirely new challenges, such as trajectory optimization for space applications.Consequently, this work aims to take the initial step in exploring QA with an entirely novel method, and the first results are definitely promising.
By leveraging the polynomial approximations, the state vector and its derivatives can be explicitly written as a linear function of the free coefficient parameter: Vector T corresponds to the i-th set of fixed polynomial coefficients, while α m = [α m x , α m y , α m z ] T represents the vector containing three independent parameters that influence the trajectory's components.By imposing the boundary conditions in Equation ( 9), we can obtain a system of m equations, where α i are considered unknowns and α m is considered a free parameter.As a consequence, all α i can be expressed as a linear function of α m and boundary conditions.By substituting the obtained expression α i in Equation ( 9), it becomes possible to obtain a linear expression as the one in Equation (10), where α m is isolated, and M(t) and N(t) are the expressions of two generic matrices as a function of time.

s(t, α
By using the general expression linear dynamical in Equation ( 8), it becomes feasible to derive control input u(t): Control vector u(t) appears only in the second term of Equation (11), where I represents the identity matrix, while Λ and H are constant matrices.These matrices are typically specific to the dynamics of the problem at hand, but in the current discussion, they are considered as generic expressions without specific values or definitions.The explicit expression of u(t) can be formulated as follows: Then, by substituting Equation (10) in Equation ( 12), control input u(t) is expressed as a linear function of α m .However, from this point on, variable α m is represented with symbol α to simplify the notation. where: where P and R are functions of variables M, N, L, and H, and their expressions depend on both the system's dynamics and the specific polynomial chosen to approximate the trajectory.As of now, they are represented in a general form, but in practical applications, these variables are given explicit expressions tailored to the specific problem at hand.Let us suppose now that our cost function aims to minimize control, and we discretize the trajectory by generating vector of d time points t = (t 0 , t 1 , . . ., t i , . . ., t d−1 , t d ).As a result, the cost function transforms into a summation of the square of the control for each time instant, and it can be expressed as follows: Value (∑ d i=0 R T i R i ) can be safely disregarded since it remains constant.Its presence in a cost function is, therefore, inconsequential, as adding or subtracting a constant value from a cost function has no impact on the result.Then, by substituting P tot = (∑ d i=0 P T i P i ) and V tot = (∑ d i=0 R T i P i ) in Equation ( 14), the cost function can be succinctly formulated as The next step to achieve the QUBO Formulation is to express each real variable in α by a number of n b binary variables.In executing this, the fixed-point binary representation is sought such that the maximum and minimum real values it can represent correspond to the boundaries of α.Indeed, referring to Equation (16), the real variable α is equal to α = α min when every z k = 0, and α = α max , every z k = 1.
where α min and α max are the upper and lower boundaries, respectively, matrix G has the size of n real × n b n real (n real = 3), g has the size of n b and represents the binary vector used for the binary form representation of integers, defined as g = (2 0 , 2 1 , . . ., 2 n b −1 ).Additionally, ∆a i = a max i − a min i (where i = x, y, z) denotes the difference between the upper and lower limits of the coefficient for each respective component.By substituting Equation (16) into Equation ( 15), an expression of the cost function in terms of binary variables z is achieved (Equation ( 18)), where L is the matrix of quadratic terms, m is the vector of linear terms, and N is a constant.By exploiting the fact that, for binary variables, z T z = z, the linear term is included into the quadratic matrix, leading to the formulation on the right hand side of Equation ( 18), with Q = Q + diag(L); furthermore, constant N can be ignored for the purpose of optimization.
The above-described methodology allows to formulate a linearized trajectory optimization as a QUBO problem, with the objective of minimizing the control employed and including the dynamical and boundaries constraints.The proposed approach offers a versatile and comprehensive solution for addressing a wide spectrum of linear optimal control problems.By transcending specific case limitations, it introduces a unified strategy that is adaptable to diverse challenges.

Lunar Landing and Rendezvous Applications
The main objective of both of these mission applications is to navigate the spacecraft in such a way that it can be safely and accurately landed or docked at a specified predetermined location.In achieving this goal, it is crucial that the spacecraft reaches this final position with specific attributes simultaneously satisfied.These attributes include having both zero relative velocity and zero relative acceleration.This ensures a controlled and precise touchdown or a docking operation, minimizing the potential for any sudden or undesirable movements during the final phases of the mission.This objective can be cast in the form of a boundary condition problem.Notably, a total of five boundary conditions (m = 5) prove sufficient to comprehensively define the entirety of this problem: Consequently, the trajectory can be approximated by a fifth-degree polynomial (as elaborated in Section 3).
In (20), five coefficients of each component are explicitly determined, while one coefficient (for each component) remains unconstrained and serves as the adjustable parameter (in this instance, α 2 = [α 2 x , α 2 y , α 2 z ]), representing the variable to be optimized through QA.By including Equation (19) in Equation ( 20), the coefficient values can be expressed as a function of boundary conditions and the free parameters: where ∆s = s f − s 0 represents the difference between the terminal and initial vector positions of the trajectory.By including Equation ( 21) into Equation ( 20), the spacecraft's trajectory can be explicitly represented as a function of α 2 . where:

Lunar Landing Problem
The mathematical expressions governing the dynamics of a lunar lander can be derived using Newton's law in the context of no central forces.In this context, the gravitational influence and the propulsive thrust produced by the propulsion system act as the only forces of motion [37].Furthermore, a flat ground is considered, thus obtaining linear dynamics.This preference is logically valid, since the present study limits its analysis to only the terminal phase of the landing trajectory.Accordingly, the formulation of the dynamics can be written as follows: where the unvarying gravitational acceleration vector of the Moon, g M = [0, 0, −g M ] T , exclusively acting along the z-axis (g M = 1.62509 m/s 2 ), along with control thrust and acceleration vectors T and u, respectively, contribute to the movement.Upon obtaining the second derivative of Equation ( 22), it becomes evident that solely time vector T undergoes differentiation, while the other terms remain held constant.The resulting s(t, α 2 ) can be included in Equation ( 23) to find the expression of control vector u: where R = N − g M represents the component that remains independent of the variable optimization.Then, by following the steps reported from Equation (15) to Equation (18) the QUBO problem can be derived.

Rendezvous Problem
The Clohessy-Wiltshire equations are a set of linearized equations commonly used in orbital mechanics to describe the relative motion between two objects in close proximity during a rendezvous or docking maneuver.These equations provide an approximation of the relative dynamics in a simplified form, making it easier to analyze and plan maneuvers [38].The Clohessy-Wiltshire equations are particularly useful when the relative distances between the two objects are small compared to the orbital distances.They are derived by linearizing the equations of motion around a reference orbit.The equations are typically used for missions where the spacecraft is performing relatively short-duration maneuvers, such as rendezvous, docking, or formation flying.The basic form of the Clohessy-Wiltshire equations is as follows: where n = µ/R is the mean motion of the reference orbit, R is the radius of the target body's circular orbit, and µ is the standard gravitational parameter.In this application, the motion is around the Earth, thus µ = 398,600 km 3 /s 2 .From Equation ( 22), the polynomial approximation and its derivatives can be included in the dynamical formulation (Equation ( 25)) to find the expression of control vector u: The control vector (in Equation ( 26)) is a linear function of variable Bα 2 and it has the same form as provided in Equation (10).Thus, by following the transposition procedure from Equation (14) to Equation ( 18), once again, the QUBO Formulation of the problem can be derived.

Results
In this section, the empirical results obtained using the D-Wave Advantage 4.1 QPU through the D-Wave cloud service "Leap" are presented.This quantum annealer adheres to a topology named Pegasus, featuring a pool of approximately 5000 qubits for computational tasks.Within this topology, each qubit establishes connections with fifteen others, resulting in a total of around 35,000 interconnections.The main objective is to assess the QA's performance in terms of computational time and optimized cost function.This evaluation is carried out through a comparative analysis with established optimization techniques.Specifically, the Quantum Annealer's performance is compared with the Interior-Point (IP) optimization method provided by the "fmincon" function in MATLAB R2022b.Additionally, a comparison is made with the genetic algorithm, Particle Swarm Optimization (PSO), which is implemented using the "particleswarm" function on MATLAB [39].Both the IP and the PSO method are implemented on a workstation equipped with 32 GB of RAM, an Intel Core i7-12700H 2.70 GHz CPU from the 12th Generation, and an NVIDIA GeForce RTX 3080 Ti GPU with 16 GB of dedicated RAM (manufactured by Intel and NVIDIA Corporation, Santa Clara, CA, USA).
The parameters used to set the simulations are provided in Table 1.We note that the variables within the problem being examined are considerably fewer in number than what the physical constraints demand.Indeed, by choosing to use 32 bits to convert the three real components into binary, the number of variables of the problem are kept very low: 3 × 32, which amounts to 96 variables, while the non-zero components of the matrix Q count 1488.Regarding the embedding procedure, it is carried out automatically by the built-in algorithms of D-Wave.The resulting embeddings are easily found and have chain lengths of eight components at maximum.As a result, the question of embedding does not arise as a concern and is resolved indirectly through the application of a substantially smaller-sized problem.The decision to limit the qubit count to 32 is not arbitrary but based on a systematic analysis of the cost functions in the context of both Lunar Landing and Rendezvous applications.Indeed, the number of qubits directly impacts the precision of conversion from real to binary numbers, as indicated in Equation (17).The outcomes of this analysis are showcased in Figure 2, where the cost functions are plotted against the base-2 logarithm of the number of qubits (log 2 (n b )).This visualization simplifies the presentation of results and underscores that a plateau is reached for both applications when log 2 (n b ) equals five, corresponding to n b = 32.As such, pursuing a greater number of variables proves to be superfluous.In addition, a Monte Carlo simulation is carried out for each test case, utilizing 10,000 samples to provide a comprehensive evaluation of QA in comparison to IP and PSO, with a focus on assessing the achieved cost function and computational time.This extensive analysis allows for a thorough understanding of the relative performance of the methods in terms of optimization quality and efficiency.

Lunar Landing Results
In Table 2, the lunar landing parameters are provided.Specifically, the initial position is selected to uphold the validity of the approximation applied to the dynamics equations, which assumes the lunar surface to be flat.As a result, the spacecraft is situated at an altitude of 3 km.This initial configuration is based on a previously published works from 2022 [36,40], ensuring a consistent and valid reference point for our analysis.The QPU results are shown in Figure 3, in terms of position, velocity, and acceleration.The figures readily demonstrate that boundary conditions are fully met; notably, the final position is reached with zero velocity and acceleration.This implies that the QA successfully yields an optimal solution for the problem at hand.The trajectory profiles and the control law are visually presented in Figures 4 and 5, featuring a comparative analysis with trajectories and control generated using MATLAB's PSO and IP methods, respectively.In both figures, the substantial overlap between trajectories offers compelling evidence of the quality of QA's solutions.Additionally, on the right side of both figures, the alignment of QA's control trajectory with those of the IP and PSO methods is illustrated, further affirming their striking similarity.This robustly demonstrates QA's effectiveness in tackling these specific problems.Furthermore, a Monte Carlo simulation is executed employing 10,000 samples to comprehensively evaluate QA's performance in contrast to IP and PSO methods.The Monte Carlo comparison results are visually represented in Figure 6, encompassing computational time and cost function optimization.In particular, the upper row reports the comparison between QA and PSO, while the bottom row reports the comparison between QA and IP methods.An insightful observation can be made when analyzing the results presented in Figure 6.Although the average value of the cost function is similar to those of QA (approximately −362.19), the PSO method (approximately −362.13), and the IP method (approximately −362.11), the QA method seems to exhibit a slightly better performance.Moreover, there is less data dispersion in QA's results.This indicates a more deterministic outcome compared to those of the other two methods.A significant difference is evident in terms of computational efficiency.Notably, the computational time required for QA (0.015 s) is an order of magnitude smaller than that of PSO (approximately 0.30 s) and IP methods (approximately 0.62 s).

Rendezvous Results
Table 3 provides the parameters for the rendezvous scenario.Notably, the final conditions are uniformly set to zero, given the utilization of equations that pertain to relative motion.It is worth noting that this scenario was taken from the study detailed in [41], where the same polynomial approximation method is used to address a rendezvous problem.This ensures consistency and aligns with the methodology established in the scientific literature.In particular, the radius of the target spacecraft orbit is set to 773 km.
The outcomes in terms of position, velocity, and acceleration are showcased in Figure 7. Similar to the previous application, the rendezvous scenario successfully achieves the targeted boundary conditions, with all values equating to zero.
Figures 8 and 9 provide a comparative analysis with trajectories and control laws generated using MATLAB's IP and PSO methods, respectively.Once again, it is evident that the results obtained from QA are highly comparable to the outcomes achieved through IP and PSO methods.Additionally, a Monte Carlo simulation is conducted using 10,000 samples to thoroughly assess QA's performance in comparison to that of IP and PSO methods.The Monte Carlo comparison results, presented in Figure 10, encompass computational time and cost function optimization.Specifically, the upper row presents the comparison between QA and PSO, while the bottom row provides the comparison between QA and the IP method.Once again, the average values of the cost function are remarkably consistent across all three methods, with QA yielding approximately −0.370, IP yielding −0.361, and PSO yielding −0.361.Furthermore, there is notably less data dispersion in the results obtained using the QA method, indicating a higher degree of determinism compared to the other two approaches.Moreover, QA demonstrates an impressively low computational time of 0.023 s, even when applied in the context of the Rendezvous scenario.This computational efficiency is notably superior, being an order of magnitude lesser than that of the IP and PSO methods, which require approximately 0.141 and 0.226 s, respectively.
Upon analyzing the results presented across different scenarios, a consistent trend emerges.While the average cost function values for QA and fmincon demonstrate comparability, a remarkable supremacy of QA in terms of computational efficiency emerges.

Conclusions
This methodology leverages QA through a QUBO translating trajectory optimization problems into the binary optimization domain.The proposed procedure presents a general approach that can be applied to optimization problems with linear dynamics.Its strength lies in its versatility, making it suitable for a wide range of scenarios within the realm of linear applications.Although it is designed for linear systems, its adaptability to different problems is a noteworthy advantage.Future research endeavors might focus on developing methods to tackle optimization challenges posed by nonlinear systems.Moreover, while the applications explored in this study may appear straightforward, a notable advantage in computational time is already discernible with the QA.This early advantage lays the groundwork for a deeper exploration of the potential harbored by this novel technology.The significance of the computational time disparity is indeed striking.The QA demonstrates a remarkable level of computational efficiency, evident through its computational time, which is orders of magnitude smaller than that of IP and PSO methods.This discrepancy in computation time may be a very important finding that could highlight an innate advantage of quantum computing in its ability to efficiently address optimization challenges.By exploiting the principles of quantum mechanics, QA is able to explore solution spaces in a parallel and probabilistic manner, enabling it to navigate through a vast number of potential solutions more quickly.In contrast, classical optimization methods often rely on iterative and sequential procedures that require greater computational efforts.
Nevertheless, it is worth noting that due to the polynomial-based approximation of the system state, an optimally exact solution might not be achieved.This approach strikes a balance between optimizing for fuel consumption and computational efficiency.Opting for higher-degree polynomials could offer more flexibility for optimization, potentially yielding more accurate outcomes.However, this decision introduces a trade-off by potentially prolonging computational time, a factor undesirable in applications requiring swift computations.Hence, the selected approach reflects a compromise between result optimality and computational efficiency.Another limitation is that the flight time, ∆t f , has to be fixed.This is, again, due to the need for reducing the optimization into a QUBO form.Indeed, if t f were a variable, it would introduce a multiplication factor with α in Equation (10).Consequently, non-linear terms would emerge, potentially reaching higher than second-degree terms within the cost function.This intricacy would render the transcription into a QUBO problem infeasible.
Given these encouraging outcomes and the numerous achievements in recent times, quantum computation has undeniably demonstrated its potential.Future research could enhance the methodology by exploring advanced trajectory approximation techniques to achieve greater precision while maintaining computational efficiency, such as B-Spline or Neural Network.In addition, the development of methods to handle nonlinear systems and to include time as a problem variable could expand the applicability of the methodology.Lastly, it might also be of interest to investigate the use of hybrid optimizers that combine quantum and classical optimization techniques, which could provide better solutions.

Figure 1 .
Figure 1.Scheme of a simple embedding procedure.

Figure 2 .
Figure 2. Lunar Landing and Rendezvous Cost Functions as a function of Number of Quibit.

Table 2 .
Parameters to set the Lunar Landing problem.

Figure 4 .
Figure 4. Comparison of trajectory and control law for Lunar Landing application between Quantum Annealing (blue line) and Particle Swarm Optimization (red line).

Figure 5 .
Figure 5.Comparison of trajectory and control law for Lunar Landing application between Quantum Annealing (blue line) and Interior Point (gree line).

Figure 6 .
Figure 6.Comparison of Monte Carlo results for Lunar Landing application between QA (blue), PSO (red), and IP methods (green) in terms of computation time and cost function.

Figure 10 .
Figure 10.Comparison of Monte Carlo results for Rendezvous application between QA (blue), PSO (red), and IP methods (green) in terms of computation time and cost function.

Table 1 .
Parameters to set the QUBO problem.

Table 3 .
Parameters to set the Rendezvous problem.