Trajectory optimization using quantum computing

TL;DR

This paper introduces a novel framework for trajectory optimization that discretizes both independent and dependent variables, enabling the use of classical and quantum search algorithms for efficient global solutions.

Contribution

It presents a new discretization scheme that reduces computational costs and integrates quantum algorithms for trajectory optimization, demonstrating potential quadratic speed-up.

Findings

Quantum search algorithms offer quadratic speed-up over classical methods.

Discretization of both variables enables flexible and efficient optimization.

Comparative analysis shows quantum methods are promising for complex trajectory problems.

Abstract

We present a framework wherein the trajectory optimization problem (or a problem involving calculus of variations) is formulated as a search problem in a discrete space. A distinctive feature of our work is the treatment of discretization of the optimization problem wherein we discretize not only independent variables (such as time) but also dependent variables. Our discretization scheme enables a reduction in computational cost through selection of coarse-grained states. It further facilitates the solution of the trajectory optimization problem via classical discrete search algorithms including deterministic and stochastic methods for obtaining a global optimum. This framework also allows us to efficiently use quantum computational algorithms for global trajectory optimization. We demonstrate that the discrete search problem can be solved by a variety of techniques including a deterministic exhaustive search in the physical space or the coefficient space, a randomized search algorithm, a quantum search algorithm or by employing a combination of randomized and quantum search algorithms depending on the nature of the problem. We illustrate our methods by solving some canonical problems in trajectory optimization. We also present a comparative study of the performances of different methods in solving our example problems. Finally, we make a case for using quantum search algorithms as they offer a quadratic speed-up in comparison to the traditional non-quantum algorithms.

Tables

Table 1. Table 1: Root mean square and the percentage errors between the minimum time obtained in Example A, B and C and the correct analytical minimum time for the brachistochrone problem.

	Root Mean Square and Percentage Error
	Example A	Example B	Example C
RMSE:	$0.0431602981808$	$0.0441659854713$	$0.00839381927831$
% Error:	$0.6$	$0.6$	$0.08$

Table 2. Table 2: Uniformly spaced sample points in x 𝑥 x for the brachistochrone problem.

	Example A
x:	$0$	$\frac{\pi }{5}$	$\frac{2\pi }{5}$	$\frac{3\pi }{5}$	$\frac{4\pi }{5}$	$\pi$
y:	$2$	$0.95$	$0.50$	$0.20$	$0.05$	$0$

Table 3. Table 3: Non-uniformly spaced sample points in x 𝑥 x for the brachistochrone problem.

	Example B
x:	$0$	$0.5$	$1.0$	$2.25$	$2.5$	$\pi$
y:	$2$	$1.05$	$0.60$	$0.10$	$0.05$	$0$

Table 4. Table 4: A finer grid resolution for obtaining a more accurate solution for the brachistochrone problem.

	Example C
x:	$0$	$0.05$	$0.25$	$0.5$	$0.75$	$1.0$	$1.5$	$2.0$	$2.25$	$2.5$	$2.75$	$2.9$	$3.05$	$3.14$
y:	$2$	$1.78$	$1.39$	$1.07$	$0.84$	$0.65$	$0.36$	$0.16$	$0.1$	$0.05$	$0.02$	$0.01$	$0$	$0$

Table 5. Table 5 : A comparison of performances of different search methods for the brachistochrone problem. (Assume ϵ = 2.46 italic-ϵ 2.46 \epsilon=2.46 .)

Search	Search method description	Cost	Minimum time	Error %
method
I	Classical exhaustive	2825761	1.0095	0.6
	(Example A, physical space)
I	Classical exhaustive	85562001	1.0107	0.7
	(Example D, coefficient space)
II	Classical randomized	5000	1.0099	0.7
III	Classical hybrid	8249	1.0085	0.5
IV	Quantum exhaustive	$\sqrt{2825761}ϵ\approxeq 4135$	1.0095	0.6
V	Quantum random	$\sqrt{5000}ϵ\approxeq 174$	1.0099	0.7
VI	Quantum hybrid	$\left(\sqrt{5000}+\sqrt{3249}\right)ϵ\approxeq 314$	1.0085	0.5

Table 6. Table 6: Uniformly spaced 4 4 4 sample points in θ 𝜃 \theta for the isoperimetric problem.

	Example G
$\theta$:	$\frac{\pi }{2}$	$\frac{5\pi }{8}$	$\frac{3\pi }{4}$	$\frac{7\pi }{8}$
$r$:	$0$	$0.26$	$0.47$	$0.62$

Table 7. Table 7: Uniformly spaced 6 6 6 sample points in θ 𝜃 \theta for the isoperimetric problem.

	Example H
$\theta$:	$\frac{\pi }{2}$	$\frac{7\pi }{12}$	$\frac{2\pi }{3}$	$\frac{3\pi }{4}$	$\frac{5\pi }{6}$	$\frac{11\pi }{12}$
$r$:	$0$	$0.17$	$0.33$	$0.47$	$0.58$	$0.64$

Table 8. Table 8: A finer grid resolution for obtaining a more accurate solution for the isoperimetric problem.

	Example I
$\theta$:	$\frac{\pi }{2}$	$\frac{9\pi }{16}$	$\frac{5\pi }{8}$	$\frac{11\pi }{16}$	$\frac{3\pi }{4}$	$\frac{13\pi }{16}$	$\frac{7\pi }{8}$	$\frac{15\pi }{16}$
$r$:	$0$	$0.13$	$0.26$	$0.37$	$0.47$	$0.55$	$0.62$	$0.65$

Table 9. Table 9 : A comparison of performances of different search methods for the isoperimetric problem. (Assume ϵ = 2.46 italic-ϵ 2.46 \epsilon=2.46 .)

Search	Search method description	Cost	Maximum area	Error %
method
I	Classical exhaustive: Example G	${101}^4$	0.174442279371647	0.052
I	Classical exhaustive: Example H	${101}^6$	0.174481616034558	0.029
I	Classical exhaustive: Example I	${101}^8$	0.174531915079274	0.0006
IV	Quantum exhaustive: Example G	${101}^2ϵ$	0.174442279371647	0.052
IV	Quantum exhaustive: Example H	${101}^3ϵ$	0.174481616034558	0.029
IV	Quantum exhaustive: Example I	${101}^4ϵ$	0.174531915079274	0.0006

Equations

$$\mathcal{J}(\vec{U}(.),\vec{X}(.),{\tau }_f)=\mathcal{M}(\vec{X}({\tau }_f),{\tau }_f)+{\int }_{{\tau }_0}^{{\tau }_f}\mathcal{L}(\vec{U}(\tau ),\vec{X}(t),t)dt.$$

$$f_l\le f(\vec{U}(\tau ),\vec{X}(t),\overrightarrow{X^{{}^{\prime }}}(t),t)\le f_{u}t\in [{\tau }_0,{\tau }_f].$$

$$h_l\le h(\vec{X}({\tau }_0),\vec{X}({\tau }_f),{\tau }_f-{\tau }_0)\le h_u,$$

$$g_l\le g(\vec{U}(t),\vec{X}(t),t)\le g_u,$$

$$\vec{X}(t)=\sum _{j=0}^{\eta -1}{\mathbf{\chi }}_j(t)$$

$${\mathbf{\chi }}_j(t)=\{\begin{array}{ll}\sum _{k=0}^{m_j}{\vec{a}}_{j,k}{\varphi }_{j,k}(t)& \text{if }t\in [t_j,t_{j+1}),\\ & \\ 0& \text{otherwise}.\end{array}$$

$$\underset{t\to t_{j+1}}{\lim }{\mathbf{\chi }}_j(t)={\mathbf{\chi }}_{j+1}(t_{j+1})\text{for }j=0,1,2\dots \eta -1.$$

$$\vec{U}(t)=\sum _{j=0}^{\eta -1}{\mathbf{\mu }}_j(t)$$

$${\mathbf{\mu }}_j(t)=\{\begin{array}{ll}\sum _{k=0}^{n_j}{\vec{b}}_{j,k}{\psi }_{j,k}(t)& \text{if }t\in [t_j,t_{j+1}),\\ & \\ 0& \text{otherwise}.\end{array}$$

$$\underset{t\to t_{j+1}}{\lim }{\mathbf{\mu }}_j(t)={\mathbf{\mu }}_{j+1}(t_{j+1})\text{for }j=0,1,2\dots \eta -1.$$

$${\text{SS}}_{j,r,k}=\{-n_{j,r,k},-(n_{j,r,k}-{ϵ}_{j,r,k}),\dots ,0,\dots ,n_{j,r,k}-{ϵ}_{j,r,k},n_{j,r,k}\}\subset \mathbb{R}.$$

$$S=\prod _{j=0}^{\eta -1}\prod _{r=1}^n\prod _{k=0}^{m_j}{S}_{j,k,r}$$

$$T=\prod _{j=0}^{\eta -1}\prod _{r=1}^m\prod _{k=0}^{n_j}{T}_{j,k,r},$$

$$O|x⟩=(-1{)}^{f(x)}|x⟩.$$

$$|\psi ⟩=H^{\otimes n}({|0⟩}^{\otimes n})=\frac{1}{\sqrt{N}}\sum _{x=0}^{N-1}|x⟩.$$

$$G=\left(H^{\otimes n}(2|0⟩⟨0|-I)H^{\otimes n}\right)O=\left(2|\psi ⟩⟨\psi |-I\right)O.$$

$$\left(2|\psi ⟩⟨\psi |-I\right)(\sum _{k=0}^{N-1}{\alpha }_k|k⟩)=\sum _{k=0}^{N-1}\left(-{\alpha }_k+2m\right)|k⟩$$

$$|\alpha ⟩=\frac{1}{\sqrt{N-M}}\sum _{x,x\text{ not marked}}|x⟩,\text{and }|\beta ⟩=\frac{1}{\sqrt{M}}\sum _{x,x\text{ marked}}|x⟩.$$

$${\tau }_f={\int }_0^{\pi }\sqrt{\frac{1+{\left(\frac{dy}{dx}\right)}^2}{2gy}}dx.$$

$${\varphi }_k(x)=\prod _{j=0,\text{leavevmode}j\ne k}^{\zeta }\frac{(x-x_j)}{x_k-x_j}.$$

$$f_A(x)$$

$$a_0$$

$$a_2$$

$$a_4$$

$$f_B(x)$$

$$a_0$$

$$a_2$$

$$a_4$$

$$a_0=2,$$

$$a_1=\frac{1}{\pi }\left(-a_0-\sum _{n=2}^{\zeta -1}-a_n{\pi }^n\right).$$

$$|y(x)|\le 2\text{for }x\in [0,\pi ].$$

$$|b_{n-2\mu }|$$

Full text

1 Introduction

The goal of trajectory optimization is to find a path or trajectory that optimizes a given quantity of interest or any other objective function associated with a certain performance measure, under a set of given constraints on the dynamics of the system. Trajectory optimization problems appear naturally in many practical situations, especially in aerospace applications. Trajectory optimization problems are important, even at a more fundamental level, as they can be related to the principle of least action or other variational principles underlying diverse physical phenomena (see feynman2005principle , goldstein2011classical , gelfand2000calculus , elsgolc2012calculus ).

Mathematically, action is a functional with the trajectory function of the system as its argument. Minimizing the action functional means selecting the optimum trajectory that minimizes the underlying integral associated to the action. Of course, in this sense, trajectory optimization is related to the subject of the calculus of variations. Therefore, the techniques proposed in this paper are also applicable to problems involving the calculus of variations.

There are many known methods of trajectory optimization (see betts1998survey , cowling2006optimal , fahroo2002direct , von1992direct ). These solution methods can be broadly classified into the categories of indirect methods and direct methods. In an indirect method one proceeds by analytically finding a set of necessary and sufficient conditions to solve a given trajectory optimization problem. Indirect methods have a long history dating back to the celebrated brachistochrone problem, posed by Bernoulli in the $16^{\text{th}}$ century (bernoulli1996problema , woodhouse1810treatise ), and which was subsequently solved by Euler, Lagrange and others by employing the techniques of calculus of variations. On the other hand direct methods have become more popular with the advent of digital computers. A direct method consists of discretizing the optimization problem and then solving the resulting non-linear optimization problem directly. Typically, in a direct method the control and the state parameters are represented by piecewise polynomial functions or a linear combination of global basis functions satisfying the boundary constraint on a discrete time grid. Several types of polynomials have been employed to represent the control and state parameters in the literature. For example in fahroo2002direct , Chebyshev polynomials are used to solve the optimization problem.

We remark that indirect methods often involve prior mathematical analysis of the problem. The resulting analytical solution is typically quite complicated and not often amenable to a computationally efficient solution. Therefore, in practical applications direct methods are often more suitable. However, in high dimensional problems direct methods also become computationally expensive. In this paper we describe a novel discretization scheme to solve the global optimization problem wherein we discretize both the independent variables (such as time) and the dependent variables. Our discretization scheme involves the selection of a finite number of coarse-grained states. Each of these states has an associated cost function. The problem then becomes a discrete search problem, in which a state with the minimum associated cost is to be determined. Then it becomes possible to employ deterministic and probabilistic methods for obtaining a global optimum. In this work, we propose a new framework for solution of the trajectory optimization problem via classical discrete search algorithms including an exhaustive search algorithm (Method I, Sect. 4.1), a random search algorithm (Method II, Sect. 4.2), and a hybrid search algorithm (Method III, Sect. 4.3). This framework also allows us to efficiently use quantum computational algorithms for global trajectory optimization. In this context, we propose new approaches for solution of the trajectory optimization problem using quantum exhaustive search algorithms (Method IV, Sect. 5.2), a quantum random search algorithm (Method V, Sect. 5.3), and a quantum hybrid algorithm (Method VI, Sect. 5.4). It turns out that quantum computers, in principle, are significantly superior to classical computers in solving the underlying discrete search problems. In fact, a slight modification of Grover’s quantum search algorithm (see grover1996fast ) leads to a solution with complexity of the order of $O(\sqrt{N})$, where $N$ depends on the number of possible states of control and state parameters. The number $N$ is more precisely defined in Sect. 3.

A main focus of this paper is to show that trajectory optimization problems can be tackled efficiently using quantum computational algorithms employed either alone or in conjunction with a randomized search. As noted earlier, to achieve this we use a discretization scheme which makes the use of quantum computing possible. We also note here that unlike many other works in the literature on trajectory optimization which are set up to only detect the local optimum, our approach enables the search of the global optimum. We will demonstrate our method using two canonical problems in trajectory optimization, namely the brachistochrone problem and the moon landing problem. The method presented here could be made even more effective if it is combined with other techniques like the gradient descent and simulated annealing (see vanderbilt1984monte ). However, those aspects are outside the scope of the present work and will be considered in a future work.

In Sect. 2, we present the problem formulation and in the following section, Sect. 3, we describe our discretization scheme. In Sect. 4 we describe classical search algorithms including the exhaustive, random and hybrid search algorithms to solve the discrete version of the trajectory optimization problem formulated in Sect. 3. Sect. 5 contains a brief introduction to quantum search methods, including Grover’s algorithm. We also describe various quantum search algorithms in this section and discuss advantages of these algorithms over their classical counterparts. Sect. 6 contains three computational examples, namely, the brachistochrone problem, the isoperimetric problem and the moon landing problem. In this section, we give several approaches to discretization and searching such as discretization in the physical space or the co-efficient space, followed by an exhaustive search, a random search, a quantum search or a hybrid search. We also discuss and compare the performances of these approaches in solving the brachistochrone problem. Finally, we present our concluding remarks in Sect. 7.

2 Trajectory optimization: Problem formulation

We will closely follow fahroo2002direct in the formulation of the trajectory optimization problem in this paper. Suppose $\vec{U}(t)\in{\mathbb{R}}^{m}$ and $\vec{X}(t)\in{\mathbb{R}}^{n}$ denote the control function and the corresponding state trajectory respectively at time $t$. The goal is to determine the optimal control function $\vec{U}(t)$ and the corresponding state trajectory $\vec{X}(t)$ for $\tau_{0}\leq t\leq\tau_{f}$ such that the following Bolza cost function is minimized:

[TABLE]

Here $\mathcal{M}$ and $\mathcal{L}$ are ${\mathbb{R}}$-valued functions. Moreover, we also assume that the system satisfies the following dynamic constraints.

[TABLE]

In addition, we also specify the following boundary conditions

[TABLE]

where $h$ an ${\mathbb{R}}^{p}$-valued function and $h_{l},h_{u}\in{\mathbb{R}}^{p}$ are constant vectors providing the lower and upper bounds of $h$. Finally, we note the mixed constraints on control and state variables

[TABLE]

with $g$ a ${\mathbb{R}}^{r}$-valued function and $g_{l},g_{u}\in{\mathbb{R}}^{r}$ are constant vectors providing the lower and upper bounds of $g$.

We note that with an appropriate transformation we may assume that $\tau_{0}=-1$ and $\tau_{f}=1$ (see fahroo2002direct ).

3 Discretization approaches

There are a number of local and global discretization methods available for the problem described in the previous section. Often the trajectory is approximated using a linear combination of a set of orthogonal polynomials. For example, one can use global spectral methods employing Chebyshev polynomials. Another possibility is to use piecewise quadratic or cubic splines to approximate the trajectory within a set of chosen discrete grid points. An appropriate choice of discretization method also depends upon the nature of the problem. For example, the Chebyshev pseudospectral method has the advantage that the node points cluster around the end points of the interval avoiding the Runge phenomenon. In this paper we present a piecewise formulation which has some advantages over global polynomial approximations. Usually, one finds in the literature that the discretization is carried out in the independent variable such as time (or space). In addition to discretizing the independent variable, in this paper we propose to discretize the dependent variable as well. In other words, discretization will mean mapping a set of continuous dependent (and independent) variables to a set of discrete dependent (and independent) variables. Next we describe our discretization scheme in detail.

We divide the time interval $[\tau_{0},\tau_{f}]$ into $\eta$ sub-intervals $[t_{j},t_{j+1}]$ for $j=0,1\cdots\eta-1$ with $\tau_{0}=t_{0}<t_{1}<t_{2}<\cdots<t_{\eta-1}<t_{\eta}=\tau_{f}$. We define

[TABLE]

with

[TABLE]

Here $\vec{a}_{j,k}\in{\mathbb{R}}^{n}$ are constant vectors and $\phi_{j,k}$ is a chosen spectral polynomial for approximating the trajectory within $t\in[t_{j},t_{j+1}]$. For example, for a quadratic approximation we may assume $m_{j}=2$ and $\phi_{j,k}=t^{k}$. Further, depending upon the nature of the problem, the following conditions may be needed to avoid any possible discontinuities

[TABLE]

Moreover, if needed, smoothness conditions involving higher-order derivatives may also be imposed.

We make a similar definition for $\vec{U}$ as,

[TABLE]

with

[TABLE]

Here $\vec{b}_{j,k}\in{\mathbb{R}}^{m}$ are constant vectors and $\psi_{j,k}$ is a chosen spectral polynomial for approximating the control function within $t\in[t_{j},t_{j+1}]$. Further, we impose the following boundary conditions similar to Eq. $\!\eqref{eq_boundary_discrete}$

[TABLE]

We note that the above formulation is more general than the one given in fahroo2002direct as on setting $m_{j}=0$ and suitably defining $\chi$ and $\mu$ we can recover the formulation given in fahroo2002direct .

Now we further discretize our problem by imposing the condition that for $1\leq r\leq n$ and for $0\leq k\leq m_{j}$ the components of vectors $\vec{a}_{j,k}(r)$ take values from a discrete set of cardinality $S_{j,r,k}$, say $\vec{a}_{j,k}(r)\in\SS_{j,r,k}$ where $\SS_{j,r,k}=\{\alpha_{1},\alpha_{2},\cdots\}$ with $\alpha_{1},\alpha_{2},\cdots\in{\mathbb{R}}$ and $\#\SS_{j,r,k}=S_{j,r,k}$. We note that the choices of $\alpha_{i}$, the elements of the set $\SS_{j,r,k}$ depend on the nature of the problem and can be appropriately modified. A typical choice could be equidistant entries of the following form

[TABLE]

Similarly, for $1\leq r\leq m$ and for $0\leq k\leq n_{j}$ we let the components of vectors $\vec{b}_{j,k}(r)$ to take values only from a discrete set of cardinality $T_{j,r,k}$, say $\vec{b}_{j,k}(r)\in{\mathbb{T}}_{j,r,k}$, where ${\mathbb{T}}_{j,r,k}=\{\beta_{i}\ |\ \beta_{i}\in{\mathbb{R}}\}$ with $\#{\mathbb{T}}_{j,r,k}=T_{j,r,k}$.

We note that with this discretization the trajectory optimization problem is now turned into a discrete search problem of size $N=ST$ where

[TABLE]

and

[TABLE]

and where the objective function to be minimized is the discrete form of Eq. $\!\eqref{Eq_Main_Cost}$. Once we take the boundary conditions given by Eq. $\!\eqref{eq_boundary_discrete}$ and Eq. $\!\eqref{eq_boundary_discrete1}$ into account, the problem size is further reduced. For example, for $j=1\ldots\eta-1$ we can always choose $\vec{a}_{j,0}$ and $\vec{b}_{j,0}$ such that the boundary conditions given by Eq. $\!\eqref{eq_boundary_discrete}$ and Eq. $\!\eqref{eq_boundary_discrete1}$ are satisfied. In any case it is clear that there are only a finite number of states, say $N$, in our search space.

In principle, we can now solve the problem by traversing our discrete search space and finding the state that corresponds to the minimum cost and satisfies all the imposed constraints. Still, in practice, the problem is the massive size of the discrete search space. We remark that even in the double precision $64$-bit representation of a variable (a real number) on a digital computer, a total of $2^{64}$ distinct states of the variable can be faithfully represented on the machine. In an optimization problem involving hundreds of state variables, clearly the size of the discrete search space could become very large and unmanageable. We propose to use a discretization scheme consisting of “coarse-grained” states, thereby reducing the size of the discrete search space compared to the traditional “fine-grained” states representation on a digital computer. The trade-off here is that reducing the size of the discrete search space affects the quality of the solution. Therefore, it becomes desirable to obtain a delicate balance between the size of the discrete search space and the quality/accuracy of the resulting solution. In most applications, the size of the discrete search space would render a direct exhaustive search computationally very costly or even outright impossible. In such cases we look for alternate practical solutions. One may use analytical techniques like the gradient method to appropriately direct the search. Another possible approach is to use a random search method to find an approximate solution and then use an exhaustive search near that approximate solution. Parallel computing may also be helpful in this context. In this article, we propose to employ a quantum search algorithm which is faster and more efficient in comparison to all the classical methods. We also discuss the possibility of using a hybrid algorithm combining random search with the quantum search. But first we consider classical approaches to solve global trajectory optimization problems.

4 New framework for trajectory optimization via classical search algorithms

In this section, we consider the problem described in Sect. 3 of finding the minimum cost, with the cost function being the discrete version of Eq. $\!\eqref{Eq_Main_Cost}$ and $N=ST$ with $S$ and $T$ as in Eq. $\!\eqref{Eq_S}$ and Eq. $\!\eqref{Eq_T}$. We will describe the application of classical algorithms, including deterministic and probabilistic approaches for solving this problem using our proposed framework of coarse-grained states. We remark that the optimum given by these methods are global, unlike many instances in the trajectory optimization literature wherein local optima are obtained.

4.1 Method I: Classical exhaustive search algorithm

We say that the search is exhaustive if each of the $N$ states have been searched for the optimum cost. An exhaustive search in the physical space is guaranteed to give an optimal global solution to any desired degree of accuracy with appropriate grid refinement. In Sect. 6.1.1 we will present further details on this method in the context of solving the brachistochrone problem.

4.2 Method II: Classical random search algorithm

In many trajectory optimization problems the discrete search space is so large that an exhaustive search becomes computationally very costly. An alternative approach in such cases is to apply the well-known “pure random search” method, wherein samples are randomly selected and the running minimum cost is updated if needed (see brooks1958discussion , patel1989pure , zabinsky1992pure ). In our context this algorithm is described below in Algorithm 1. In Sect. 6.1.3, we will present a computational example employing this algorithm for solving the brachistochrone problem.

We note that this random search algorithm is guaranteed to succeed, with probability $1$, as the number of steps $n$ approaches infinity. To make this more precise, suppose that all the states are equally likely to yield the minimum cost. Let $p(r,n)$ denote the probability of the random search algorithm picking the $r^{\text{th}}$ ranked minimum cost path in the list, in $n$ random searches. Clearly $p(r,n)=\sum_{k=0}^{n-1}(\frac{r}{N})(\frac{N-r}{N})^{k}=1-(1-\frac{r}{N})^{n}$ and $\lim\limits_{n\to\infty}p(r,n)\to 1$.

4.3 Method III: Classical hybrid search algorithm

The hybrid search algorithm is a combination of the random and the exhaustive search algorithms and in many situations it is the most efficient classical search method. In a hybrid search algorithm first a random search is carried out to obtain an approximate solution. Then the result of the random search is further improved by performing an exhaustive search in a finer discrete grid near the approximate solution given by the random search. A computational example will be given in Sect. 6.1.4.

5 New framework for trajectory optimization via quantum search algorithms

Quantum computing is an exciting field which beautifully combines quantum physics with computer science. There has been impressive progress in theoretical development of quantum algorithms as well as in their practical implementation on quantum computers. We refer the readers to any standard book on the subject for a more in-depth treatment (for eg., nielsen2002quantum , rieffel2011quantum or yanofsky2008quantum ).

5.1 Overview of quantum computation

A fundamental concept in quantum computation is that of the quantum bit, or ‘qubit’. Unlike the classical bit, the qubit, say $\ket{q}$, can be in a superposition of the states $\ket{0}$ and $\ket{1}$, i.e., $\ket{q}=a\ket{0}+b\ket{1}$ with $a,b\in{\mathbb{C}}$. Hence, a qubit can be considered a vector in a two-dimensional complex vector space with the set $\{\ket{0},\ket{1}\}$ forming a basis. A quantum system can be transformed from one state to another by the so called unitary transformations, which are reversible (see Fig. 1, for some examples). Once a classical measurement is taken the state of the qubit $\ket{q}$ is changed from the superposition of the states $\ket{0}$ and $\ket{1}$ to either [math] or $1$, with respective probabilities of $|a|^{2}$ and $|b|^{2}$. This measurement step is a non-unitary and non-reversible transformation.

We give a few examples to illustrate the working of a single qubit system. So, mathematically, if the state of a single qubit $\ket{q}=a\ket{0}+b\ket{1}$ is represented by the matrix $\left[\begin{smallmatrix}a\\ b\end{smallmatrix}\right]$ and if a unitary transformation is represented by the $2\times 2$ unitary matrix $U=\left[\begin{matrix}[r]\alpha&\beta\\ \gamma&\delta\\ \end{matrix}\right]$ then the resulting state is given by $Uq=\left[\begin{matrix}[r]\alpha&\beta\\ \gamma&\delta\\ \end{matrix}\right]\left[\begin{smallmatrix}a\\ b\end{smallmatrix}\right]$. Suppose the unitary matrix is $X=\left[\begin{matrix}[r]0&1\\ 1&0\\ \end{matrix}\right]$ then the result of action of the unitary transformation $X$ is $\left[\begin{matrix}[r]0&1\\ 1&0\\ \end{matrix}\right]\left[\begin{smallmatrix}a\\ b\end{smallmatrix}\right]=\left[\begin{smallmatrix}b\\ a\end{smallmatrix}\right]$. The unitary matrix $X$ represents the quantum NOT gate as it swaps the computational basis, i.e., it sends $\ket{0}\to\ket{1}$ and $\ket{1}\to\ket{0}$, see Fig. 1. Another very useful quantum gate is the Hadamard gate, whose matrix is given by $H=\frac{1}{\sqrt{2}}\left[\begin{matrix}[r]1&1\\ 1&-1\\ \end{matrix}\right]$. It is often used to condition the input qubit which is often $\ket{0}$ to a uniform superimposed state, $H\ket{0}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$. The above discussion can be easily extended to the multiple qubit systems. For example, a two qubit system has $4$ possible states represented by $\ket{00},\ket{01},\ket{10}$ and $\ket{11}$.

5.1.1 **Quantum search: Grover’s algorithm **

We will briefly describe Grover’s quantum search algorithm grover1996fast for the sake of clarity and continuity of presentation in this work. Suppose that there is a list of $N$ items, and we are required to search for one or more marked items. Also assume that we have access to an oracle, who when presented with an item from the list will answer if the item is marked or not. For example, the oracle can be assumed to be a black box function such that for any item with index $x$ in the list, i.e., with $0\leq x\leq N-1$, we have $f(x)=1$ if $x$ is a marked item, otherwise $f(x)=0$. It is clear that with the classical computation total $O(N)$ calls to oracle will be required for the search problem. However, with the quantum computation only $O(\sqrt{N})$ oracle calls are needed to solve the search problem with probability $1$. Furthermore, if there are $M$ solutions of the search problem (of size $N$) then only $O(\sqrt{N/M})$ calls to the oracle is sufficient to solve the search problem.

For simplicity, let us assume that $N=2^{n}$. Therefore, the discrete search space can be represented by a $n$ qubit system. The algorithm starts with the system in $\ket{0}^{\otimes n}$ state, i.e., all the $n$ qubits are initially in the state $\ket{0}$. The action of the oracle $O$ (see Eq. $\!6.3$ in nielsen2002quantum ) on the input $\ket{x}$ can be represented by

[TABLE]

The action of the Hadamard transform $H^{\otimes n}$ on the input $\ket{0}^{\otimes n}$ is used to transform the input to a uniform superimposed state, say $\ket{\psi}$, such that

[TABLE]

Next the algorithm proceeds with the repeated application of Grover operator $G$ which is shown in the Fig. 2. The third step in the Fig. 2 is the phase shift which is transforming all non-zero basis state $\ket{x}$ to $-\ket{x}$ with $\ket{0}$ remaining unchanged. It is easy to see that the phase shift $\phi$ is the unitary operator $(2\ket{0}\bra{0}-I)$. It follows that the Grover operator $G$ is essentially the unitary transformation

[TABLE]

We note that at the heart of Grover’s algorithm lies the operator $G$. In fact, as noted earlier $G$ consists of $\left(2\ket{\psi}\bra{\psi}-I\right)O$. The action of oracle amounts to changing the phase of the marked item. On the other hand it can be checked that

[TABLE]

where $m=\frac{1}{N}\sum_{k=0}^{N-1}\alpha_{k}$ is the mean of $\alpha_{k}$. Hence, it is clear that the action of $\left(2\ket{\psi}\bra{\psi}-I\right)$ is essentially an inversion about the mean. Therefore, one iteration of Grover operator results in amplification of the amplitude of the marked item (See Fig. 3).

We note that if $M=1$ then the Grover operator $G$ is applied $R\approx\frac{\pi}{4}\sqrt{N}$ times. After $R$ iteration of $G$ the measurement of the system yields the marked item. To summarize Grover’s algorithm for the case $M=1$ is as follows.

Initialize the input $n$ qubits as $\ket{0}^{\otimes n}$. 2. 2.

Apply $H^{\otimes n}$ 3. 3.

Apply Grover operator $G$ a total of $R\approx\frac{\pi}{4}\sqrt{N}$ times 4. 4.

Measure the $n$ qubits.

Some remarks:

It is important to note that there is an optimum number of times (the integer $R$ as above) that the Grover’s iteration should be carried out to obtain the maximum probability of success (close to $1$). This is in contrast to the classical behavior where the more iterations lead to better results. 2. 2.

The Grover operator $G$ (see Eq. $\!\eqref{Eq_Grover_operator}$) could be very nicely explained as rotation in a two dimensional space spanned by the vectors $\ket{\alpha}$ and $\ket{\beta}$ with

[TABLE]

The initial state $\ket{\psi}$ is in the span of $\ket{\alpha}$ and $\ket{\beta}$, as it is easily verified that $\ket{\psi}=\sqrt{\frac{N-M}{N}}\ket{\alpha}+\sqrt{\frac{M}{N}}\ket{\beta}$. If we set $\cos(\frac{\theta}{2})=\sqrt{\frac{N-M}{N}}$ then the Grover operator may be written as $\left[\begin{matrix}[r]\cos(\theta)&-\sin(\theta)\\ \sin(\theta)&\cos(\theta)\\ \end{matrix}\right].$ The number of iteration $R$ is chosen such that $R\theta$ is the angle required to turn the state vector $\psi$ very close to $\ket{\beta}$. For example, if $1=M<<N$ then we get $\theta\approx\sin(\theta)=2\sin(\frac{\theta}{2})\cos(\frac{\theta}{2})\approx 2\sqrt{\frac{1}{N}}$. Also as $\psi$ is initially almost along the $\ket{\alpha}$ the total turn is approximately $\frac{\pi}{2}$. Therefore, we must have $R\theta\approx 2R\sqrt{\frac{1}{N}}\approx\frac{\pi}{2}$. And we see that $R\approx\frac{\pi}{4}\sqrt{N}$. We refer readers to nielsen2002quantum for further details.

5.2 Method IV: Quantum exhaustive search algorithm

Several quantum search algorithms based on Grover’s algorithm have been treated in the literature. If one is looking for $M$ marked things out of $N$ things, with number $M$ not known a priori, then the expected running time of these algorithms are $O(\sqrt{\frac{N}{M}})$. Dürr and Høyer have given a quantum algorithm in durr1996quantum to find the minimum cost. Algorithm 2 described below is an adaptation of their algorithm suited for our framework to solve a trajectory optimization problem. We will give a computational example in Sect. 6.1.5. We note that Algorithm 2 finds the global minimum for the problem described in Sect. 3 with the cost function being the discrete version of Eq. $\!\eqref{Eq_Main_Cost}$. It takes a list of size $N$ as input. This list contains all the possible states of the system. The algorithm proceeds by first selecting a random state from the input list and setting its associated cost as the running minimum. The value of parameter $\lambda$ and the variable $m$ is also initialized. The variable $m$ is initially set to $1$ and later in the algorithm it is adaptively scaled by $\lambda$, by setting its value to $\lambda m$ if the search for a state with lower associated cost than the current running minimum cost is not successful. The variable $m$ effectively controls the number of Grover’s rotation $r$ that is performed. A priori, it is not known that how many of the states have their associated costs less than the current running minimum cost and therefore $m$ is adaptively scaled to ensure that a correct value for $r$ is selected as required by Grover’s search algorithm.

The algorithm runs until the total accumulated number of Grover’s rotation $G$ is less than $22.5\sqrt{N}+1.4\log^{2}(N)$. We refer readers to durr1996quantum for further details. We also note that the running time of the above algorithm is optimal up to a constant prefactor. See bulger2003implementing and lara2014new for other variants of the above algorithm with smaller constant prefactors. We note that in place of Algorithm 2, these other variants may also be employed to solve the trajectory optimization problem.

It is clear that the quantum algorithms discussed above are of the order of $O(\sqrt{N})$. Hence, they are far superior to their classical counterparts which are of the order of $O(N)$.

5.3 Method V: Quantum random search algorithm

As before we consider the problem of finding the minimum cost for the problem described in Sect. 3, with the cost function being the discrete version of Eq. $\!\eqref{Eq_Main_Cost}$ and $N=ST$ with $S$ and $T$ as in Eq. $\!\eqref{Eq_S}$ and Eq. $\!\eqref{Eq_T}$. As the first step of this algorithm a predetermined number of states are selected uniformly at random out of the total $N$ states. Then Algorithm 2 is applied on these selected states. Similar to the classical random search algorithm, this algorithm is especially useful in situations wherein the discrete search space is so big that an exhaustive search has a prohibitively high computational cost. We will consider a computational example in Sect. 6.1.6.

5.4 Method VI: Quantum hybrid search algorithm

We consider a quantum hybrid algorithm to solve the problem described in Sect. 3 with the cost function being the discrete version of Eq. $\!\eqref{Eq_Main_Cost}$ and $N=ST$ with $S$ and $T$ as in Eq. $\!\eqref{Eq_S}$ and Eq. $\!\eqref{Eq_T}$. A quantum hybrid search algorithm is a combination of a quantum random search algorithm and a pure exhaustive quantum search algorithm. It combines the best attributes of both a quantum random search algorithm and a pure quantum exhaustive search algorithm. To begin with, in a hybrid algorithm first a fixed number of states are selected uniformly at random out of the total $N$ states and on these selected states the quantum algorithm described earlier (Algorithm 2) is applied. Next, using the approximate solution thus obtained a finer discretization is carried out. Finally, an exhaustive quantum search is carried out using Algorithm 2 in this finer grid. A computational example is given in Sect. 6.1.7.

Here it is interesting to note that bulger2007combining also describes a hybrid method, by effectively combining local search with Grover’s algorithm, to obtain the global optimum of an objective function. The quantum basin hopper algorithm, proposed in bulger2007combining , considers a black-box real-valued objective function $f$ defined on a discrete domain $S$, of size $2^{n}$. It also makes certain regularity assumptions on the objective function and its domain. As the method in bulger2007combining does not deal with the problem of discretization, it is not directly applicable for solving a continuous trajectory optimization problem. It is in this context that our framework could be useful. Following our discretization approach, depending upon the presence of any regularity or local structure in the problem, various local search methods and variants of Grover’s algorithms, including the method in bulger2007combining , can be used in our framework. Although, in our present work (including the computational examples that follow), we do not make any assumptions about the availability of local search methods or about the presence of any local structure in the optimization problems.

6 Computational examples

Now we consider a couple of examples for illustrating our method.

6.1 Brachistochrone problem

The first example that we take is the well-known Brachistochrone problem. The problem is to find the required trajectory of a particle starting from the rest under the influence of gravity, without any friction, such that it slides from the one fixed point to the other in the least possible time. We let $\vec{X}(t)=x(t)\vec{i}+y(t)\vec{j}$ to represent the position of the particle. The boundary condition that we consider is $(x(\tau_{0}),y(\tau_{0}))=(0,2)$ and $(x(\tau_{f}),y(\tau_{f}))=(\pi,0)$ with $\tau_{0}=0$. The goal is to minimize $\tau_{f}$ given by

[TABLE]

Next we describe some possible approaches for discretization that could be employed to solve the brachistochrone problem.

6.1.1 Global dicretization in physical space

Out of the infinitely many possible paths the goal is to pick the path with the minimum time. For practical purposes it is sufficient to find a ‘good enough’ approximate solution (based on acceptable levels of errors). The discretization is chosen based on what constitute a ‘good enough’ solution in a given context. In our present example, the physical space consists of the rectangle $[0,\pi]\times[0,2]$ in ${\mathbb{R}}^{2}$. Of course, one may as well consider a bigger rectangle to allow for the possibility of better solutions. One can now discretize the rectangle $[0,\pi]\times[0,2]$ in several possible ways. For example, let $x_{k}=\frac{k\pi}{\zeta}$ for $k=0,1\cdots\zeta$. Let $y_{k}\in\{y_{[i]}=\frac{2i}{L}:\>i=0,1\cdots L\}$ for $k=1$ to $\zeta-1$. We set $y_{0}=2$ and $y_{\zeta}=0$ to take into account the boundary condition. Next we use Lagrange interpolation to set $y(x)=\sum_{k=0}^{\zeta}y_{k}\phi_{k}(x)$. Here $\phi_{k}(x)$ is chosen such that $\phi_{k}(x_{j})=1$ if $j=k$ and $\phi_{k}(x_{j})=0$ otherwise. More explicitly,

[TABLE]

Essentially, it means that once a discretization of $x$ is carried out then we require that at $x=x_{k}$, the corresponding $y(x_{k})$ can only have values from the set $\{y_{k}\}$ (see Fig. 5). We note that with the above discretization there are $N=(L+1)^{\zeta-1}$ total possible cases need to be considered.

Theoretically, one can always approach the correct solution to the desired accuracy by making the discrete grid finer. However, the discrete search space grows very rapidly and other methods should be combined with exhaustive search in practical applications. We consider a few examples to compare the correctness of the results with respect to the size of the discrete search spaces and the points chosen for discretization. Examples A and B are chosen to highlight the effect of uniform and non-uniform grid spacing respectively (for $\zeta=5$). Example C is chosen to highlight the effect of having a finer grid resolution.

• •

Example A: Consider the path obtained by using the Lagrange interpolation on the $(x,y)$ coordinates listed in Table 2,

[TABLE]

with

[TABLE]

This path can be discovered by taking $L=40$ and $\zeta=5$ in the preceding discussion. Therefore, the total size of the discrete search space is $41^{4}=2825761$. We note that the time obtained in this case, $1.00946330885$ sec, differs from the correct analytical solutions $1.0035449615773016$ sec by only about $0.6$%.

• •

Example B: Let us now consider a different discretization and pick the points as shown in Table 3. The path obtained by using the Lagrange interpolation in this case is

[TABLE]

with

[TABLE]

This path can be discovered by using a discrete search space of size $41^{4}=2825761$. The minimum time obtained is $1.00960568357$ sec. The error in the time obtained in this case is $0.6$%, the same as the previous case.

• •

Example C: In this example we get a more accurate result by considering more points as shown in Table 4. The size of the discrete search space to discover this path, if we search exhaustively, is $201^{12}$. The minimum time obtained in this case is $1.00436810144$ sec and the error percentage is about $0.08$%.

In Fig. 6 the paths for the above examples are plotted along with the cycloid resulting from the correct analytical solution. Errors associated with these paths are shown in Table 1. We note that results shown in Fig. 6 for Examples A, B and C are in good agreement with the analytical solution. As expected, Example C has the least error (as noted in Table 1) which can be attributed to a finer grid resolution.

Some remarks:

One can easily put the above description in the general discretization framework that we discussed in the Sect. 3. For example, setting $x(t)=t$ and $y(t)=f(t)$ is the same as saying $y=f(x)$. 2. 2.

In the above discussion, we have picked the discrete search spaces rather naively, just to show that it is possible to get a good solution (within $0.6\%$ accuracy) even in such cases. One can get a smaller discrete search space by starting with an initial guess, for example in brachistochrone problem this could be a linear path, or one can do multiple passes of search by progressively refining the discrete grid in later passes.

6.1.2 Discretization in coefficient space

Sometimes due to the nature of the optimization problem it is more convenient and computationally efficient to discretize in the coefficient space. We will solve the brachistochrone problem by employing this technique. Suppose $y(x)=\sum_{n=0}^{\zeta-1}a_{n}x^{n}$ is a possible path. We want to discretize the coefficients $a_{n}$. In order to perform this discretization, we need to first determine the range of the coefficients $a_{n}$, i.e., their minimum and maximum values. The boundary condition, $y(0)=2$ implies

[TABLE]

and we use the boundary condition $y(\pi)=0$ to fix

[TABLE]

Now guided by the physical space of this problem we impose the condition that

[TABLE]

Further, we note some interesting results in approximation theory (ref. rahman2002analytic , Theorems $16.3.1$ and $16.3.2$ therein) on the bounds on coefficients of bounded polynomials. According to these theorems, given a bound on $y(x)$, such as in Eq. \!\eqref{eq_y_constraint}, the bound on the coefficients of the polynomial $y(x)$ can be determined by using the coefficients of Chebyshev polynomials of the first kind. In the following, we provide more details on these theorems to demonstrate their application in our context.

Let $T_{n}(x)$ be the Chebyshev polynomial of the first kind of degree $n$ with the series expansion given as $T_{n}(x)=\sum_{k=0}^{n}t_{n,k}x^{k}$. Let $\frac{1}{2}y(\frac{\pi}{2}(x+1))=\sum_{n=0}^{\zeta-1}b_{n}x^{n}$. It is clear that $|\frac{1}{2}y(\frac{\pi}{2}(x+1))|\leq 1$ for $|x|\leq 1$, and therefore it satisfies the hypothesis of the earlier mentioned theorems in rahman2002analytic . Hence, we obtain the following bounds on the coefficients of $\frac{1}{2}y(\frac{\pi}{2}(x+1))$

[TABLE]

Of course, the coefficients $a_{1},\cdots a_{\zeta-1}$ can be expressed in terms of the coefficients $b_{1},\cdots b_{\zeta-1}$ using the relation $\sum_{n=0}^{\zeta-1}b_{n}x^{n}=\frac{1}{2}\sum_{n=0}^{\zeta-1}a_{n}\left(\frac{\pi}{2}(x+1)\right)^{n}$, and therefore from the above relations using the bounds on coefficients $b_{1},\cdots b_{\zeta-1}$ coefficients $a_{1},\cdots a_{\zeta-1}$ can be discretized.

Example D: We set $\zeta=6$ and note that $T_{5}(x)=16x^{5}-20x^{3}+5x$ and $T_{4}(x)=8x^{4}-8x^{2}+1$. On using Eq. \!\eqref{eq_boundary_a_0} through Eq. $\!\eqref{eq_coefficient_bound2}$, we obtain the following conditions

[TABLE]

Using the above relations, we discretized the coefficients $a_{2},a_{3},a_{4}$ and $a_{5}$ by letting each of $b_{2}$ and $b_{4}$ to take values in the set $\{-8,-7.8,\cdots 7.8,8\}$ of cardinality $81$, $b_{3}$ to take values in the set $\{-20,-19.5,\cdots,\allowbreak 19.5,20\}$ of cardinality $81$, and $b_{5}$ to take values in the set $\{-16,-15.8,\cdots,15.8,16\}$ of cardinality $161$. Then one can find the path $y(x)=-0.0418x^{5}+0.3942x^{4}-1.4449x^{3}+2.7559x^{2}-3.1831x+2$ with the corresponding minimum time of $1.01067460259$ sec. We note that the error in the minimum time calculated is about $0.7\%$ and the total number of paths considered is $81\times 81\times 81\times 161=85562001$. Of course, the accuracy can be further improved by making the discretization grid finer as well as increasing $\zeta$.

A remark: In the above example, an error of $0.7\%$ in the minimum time was obtained with a search space of the size $85562001$. A much smaller search space will suffice, if a larger error in calculating the minimum time is allowed. For comparison we note that, using a similar approach as discussed above for a degree three polynomial, i.e., for $\zeta=4$, we obtained the minimum time with an error of $1.4\%$, by considering a search space of the size $1600$.

6.1.3 Classical random search algorithm in physical space

Example E: Now we consider the same global discretization using the actual physical space that we considered in the Example A, Sect. 6.1.1, with $\zeta=5$ and $L=40$. Then the best path that we obtained on a trial run was

[TABLE]

with

[TABLE]

and the corresponding minimum time found was $1.00993422472$ sec. We note that the error in the minimum time calculated is about $0.7\%$ and the total number of paths considered is $5000$. Of course, the accuracy can be increased by making the discretization grid finer as well as increasing the number of paths considered.

6.1.4 Classical hybrid search algorithm

In many situations, it is a good idea to follow the chosen discretization method by a random search to obtain an approximate solution. Then one can pick a finer discrete grid to further improve the result of the random search by performing an exhaustive search in a finer discrete grid near the approximate solution given by the random search.

Example F: The result obtained from the random search in the Sect. 6.1.3 is further refined by performing an exhaustive search near the best path given by the random search algorithm. We consider $x_{k}=\frac{k\pi}{\zeta}$ for $k=0$ to $\zeta$ and $\zeta=5$. Further, we set $y(x_{0})=2$, $y(x_{1})\in\{0.8,0.85,0.90$, $0.95\}$, $y(x_{2})\in\{0.45,0.50,0.55,0.60\}$, $y(x_{3})\in\{0.20,0.21,0.22,\cdots 0.39\}$ and $y(x_{4})\in\{0.00,0.01,0.02,\cdots 0.19\}$ and $y(x_{5})=\pi$. The best path obtained is

[TABLE]

with

[TABLE]

and the corresponding minimum time $1.00852712176$ sec. We note that the error in the minimum time calculated is about $0.5\%$ and the total number of paths considered is $5000$ for the random search and subsequently $3249$ for the exhaustive search, i.e., a total of $8249$ paths.

6.1.5 Quantum exhaustive search algorithm

As noted earlier, quantum search algorithms give a quadratic speed-up compared to the classical algorithms in search problems. Therefore, the quantum search algorithm, Algorithm 2, could be employed to perform an exhaustive search. We remark that the computational cost of Algorithm 2 in performing an exhaustive search on the states resulting from the same discretization as considered in Example A, is $4135$ whereas the cost for the equivalent classical exhaustive search was $2825761$ as noted earlier (see the Table 5).

6.1.6 Quantum random search algorithm

The random search that we considered earlier using the classical algorithm in Example E, Sect. 6.1.3, can be carried out using the quantum search algorithm. If we simply pick $5000$ random paths and apply the quantum search algorithm, Algorithm 2, on these paths to select the optimum path, then the cost of the quantum algorithm would be $\sqrt{5000}\,\epsilon\approxeq 174$ paths, where according to Theorem 6.2, bulger2003implementing , the constant $\epsilon$ has be taken to be approximately equal to $2.46$

6.1.7 Quantum hybrid search algorithm

The hybrid quantum algorithm considered here is a two step quantum algorithm. The first step of this algorithm is an application of the quantum random search, in which the quantum algorithm, Algorithm 2, is applied on a predetermined number of states selected uniformly at random out of the total $N$ states. The first step results in an approximate solution of the trajectory optimization problem. In the second step, we propose to employ Algorithm 2 again, following the random search, for performing an exhaustive search in a finer discrete grid near the approximate solution obtained earlier. Clearly such a method combines the best features of both the randomized algorithm and the exhaustive quantum search algorithm and can be advantageous in many situations. If we continue with the example of the random quantum algorithm discussed in the previous subsection as the first step of the hybrid quantum algorithm then the corresponding cost of the first step would be $\sqrt{5000}\,\epsilon\approxeq 174$. The computational cost for the subsequent step (with application of Algorithm 2 again but on the paths considered in Example F, Sect. 6.1.4) would be $\sqrt{3249}\,\epsilon\approxeq 140$. Hence, the total computational cost for the hybrid quantum algorithm is $\approxeq 174+140=314$. A summary of the costs associated to various methods is given in the Table 5.

6.2 Isoperimetric problem

As our second example, we consider the well-known isoperimetric problem bliss1925calculus . The objective in this problem is to find the maximum area enclosed between a curve of the given fixed length and a given fixed straight line ${\mathfrak{l}}$, such that the endpoints of the curve lie on the straight line ${\mathfrak{l}}$. This problem could be solved by employing the principles of calculus of variations. Here we will obtain a numerical solution of this problem based on our proposed framework.

Let $r=r(\theta)$, for $\theta=\frac{\pi}{2}$ to $\theta={\pi}$, be a representation of the curve in polar coordinates with the boundary condition $r(\frac{\pi}{2})=0$. Let the length of the curve be $\frac{\pi}{3}$, i.e.,

[TABLE]

Let the fixed straight line be the line $\theta=\pi$. The objective function to be maximized is the area $A$ given by

[TABLE]

Next we explain our discretization scheme for this problem. The variable $\theta$ is treated as the independent variable and the variable $r=r(\theta)$ is assumed to be the dependent variable. The parameter space consisting of the rectangle $[\pi/2,\pi]\times[0,b]$ in ${\mathbb{R}}^{2}$ is discretized, with a suitably chosen $b\leq\frac{\pi}{3}$. We let $\theta_{k}=\frac{k\pi}{\zeta}$ for $k=0,1\cdots\zeta$. Let $r_{k}=r(\theta_{k})\in\{r_{[i]}=\frac{bi}{L}:\>i=0,1\cdots L\}$ for $k=1$ to $\zeta-1$. We set $r(\pi/2)=0$ to take into account the boundary condition. Further, using Lagrange interpolation we can determine the function $r=r(\theta)$ as $r=c\sum_{k=0}^{\zeta}r_{k}\phi_{k}(\theta)$ with $c$ chosen such that Eq. $\!\eqref{length_condition}$ is satisfied. Here $\phi_{k}(\theta)$ is chosen such that $\phi_{k}(\theta_{j})=1$ if $j=k$ and $\phi_{k}(\theta_{j})=0$ otherwise. More explicitly

[TABLE]

Essentially, it means that once a discretization of $\theta$ is carried out then we require that at $\theta=\theta_{k}$, the corresponding $r(\theta_{k})$ can only have values from the set $\{r_{[i]}=\frac{bi}{L}:\>i=0,1\cdots L\}$. We note that with the above discretization there are $N=(L+1)^{\zeta}$ total possible paths. We fix $b=1$ and $L=100$ for the following examples.

6.2.1 Global discretization in physical space

Examples G, H and I are given below to show the effect of increasing grid resolution on the accuracy of solutions.

• •

Example G: We consider the path obtained by using the Lagrange interpolation using the data of polar-coordinates $(r,\theta)$ contained in Table 6,

[TABLE]

with

[TABLE]

This path can be discovered by considering $\zeta=3$ in the preceding discussion. The size of the discrete search space is $101^{4}$. We note that the maximum area obtained in this case, $0.174442279371647$ sq. units, differs from the correct analytical solutions $0.17453292519$ sq. units by only about $0.052$%.

• •

Example H: Let us now consider a different discretization and pick the points as shown in Table 7. The path obtained by using the Lagrange interpolation in this case is

[TABLE]

with

[TABLE]

This path can be discovered by using a discrete search space of size $101^{6}$. The maximum area obtained is $0.174481616034558$ sq. units. The error in the time obtained in this case is $0.029$%, which is better than the previous case.

• •

Example I: In this example, we get a more accurate result by considering even more points as shown in Table 8. The size of the discrete search space to discover this path, if we search exhaustively, is $101^{8}$. The maximum area obtained in this case is $0.174531915079274$ sq. units and the error percentage is less than $0.0006$%.

6.2.2 Quantum exhaustive search algorithm

The quantum search algorithm, Algorithm 2, could be employed to perform an exhaustive search using the same discretization as considered in Examples G, H and I to gain significant improvement in the size of the search space as discussed earlier. In fact the size of the search space on using Algorithm 2, turns out to be $101^{2}\,\epsilon$, $101^{3}\,\epsilon$ and $101^{4}\,\epsilon$, for the Examples G, H and I respectively, with $\epsilon=2.46$ (see the Table 9).

6.3 Moon landing problem

Our third example is the Moon landing problem considered in fahroo2002direct and murecsan2012soft . The problem is to find the optimum control for soft landing a spacecraft on the surface of the moon such that the fuel consumption is the minimum. We assume that the motion is vertical and that the lunar gravitation is a constant, $g=1.63$, throughout the motion. Let $m(t)$ be the mass of the spacecraft including the fuel at the time $t$. Let $h(t)$ and $v(t)$ denote the height and the velocity of the spacecraft at the time $t$. The thrust at time $t$, say $T(t)$, is given by $-k\frac{dm(t)}{dt}=-ku(t)$ and we treat $u(t)=\frac{dm(t)}{dt}$ as a control variable. The control problem can now be described as

[TABLE]

subject to the equations of motion

[TABLE]

and the constraints

[TABLE]

where $\mu$ is a constant which determines the maximum thrust. Also, suppose the total time of flight is $\tau$ then the condition for the soft landing is that $v(\tau)=0$ and $h(\tau)=0$ and $h(t)>0$ for $t<\tau$. The initial condition is given by

[TABLE]

It easily follows on integrating the second equation of motion in Eq. $\!\eqref{eq_motion}$ that

[TABLE]

Therefore, instead of maximizing $m(t)$ one can consider the equivalent objective function of minimizing the time of flight $\tau$.

Next we describe our method for solving this equivalent control problem. We consider the following discretization. We pick $0=t_{0}<t_{1}<t_{2}<t_{3}<t_{4}=t_{\text{max}}$ as time points for discretization of the problem. We set

[TABLE]

Here, $a_{2},a_{3}$ and $a_{4}$ are chosen from a discrete set $\{\alpha_{k}\}$ for $k=1$ to $n$ with $-\mu\leq\alpha_{k}<0$ and $a_{1}$ is to be determined later. We note that, now with the help of equation of motions and the initial conditions $m(t),v(t)$ and $h(t)$ could be completely determined. In fact we get,

[TABLE]

and

[TABLE]

Next, on integrating $v(t)$ we get $h(t)=h_{0}+\int_{0}^{t}v(x)\,dx$. We treat $a_{1}$ and $\tau$ as a variable to solve the constraint equations for a soft landing, i.e., we solve for $\tau$ and $a_{1}$ such that $h(\tau)=0$ and $v(\tau)=0$ with $h(t)>0$ for $t<\tau$. We also let $a_{2},a_{3}$ and $a_{4}$ take value in the discrete set $\{-15,-14.99,-14.98,\cdots 0.1\}$. For the initial conditions $h_{0}=50002.65,v_{0}=-178,m_{0}=2500$ and the sampling points $t_{1}=25,t_{2}=75,t_{3}=200$ and $t_{4}=400$ with $g=1.63,k=585,\mu=-15,m_{a}=800$ we obtain the following solution: $a_{1}=-0.0018,a_{2}=-6.30,a_{3}=-6.30,a_{4}=-6.35$ and $\tau=290.98$ and $m(\tau)=819.76$. The time evolution of the mass, velocity and height as well as the variation of height versus velocity are shown in figures 11 through 13. These results appear to be in agreement with the known solutions of this problem. We also note that these results are based on Method I (see Sect. 4.1). All the other methods (Method II to Method VI) may also be applied to solve this problem giving similar results. The performance comparisons of these methods are also expected to be similar those given in Table 5 for the brachistochrone problem, with quantum algorithms outperforming their classical counterparts.

Some remarks:

We have arbitrarily chosen $a_{1}$ as the variable chosen for satisfying the constraints for a soft landing. One can easily make an alternative choice. 2. 2.

The choices of the time instants $t_{1},t_{2},t_{3}$ and $t_{4}$ were also arbitrary. In fact, one can let $t_{1},t_{2},\cdots t_{4}$ take values in an appropriately chosen discrete set and then pick the best possible solutions.

7 Conclusion

We have shown that a trajectory optimization problem (or a problem involving calculus of variations) can be formulated in the form of a search problem in a discrete space. An important feature of this work is our formulation of the discretization of the optimization problem incorporating the treatment of the dependent variable space in terms of appropriately chosen coarse grained states. The number of the coarse grained states is chosen to be significantly less than the traditional digital computer representations. This not only reduces the computational cost but also enables us to employ deterministic and stochastic approaches for obtaining global optimum. In particular, we presented the use of our proposed discretization approach and classical methods (Methods I–III) to solve the trajectory optimization problem. The proposed framework also enables the use of quantum computational algorithms for global (trajectory) optimization. In this work, we showed that the discrete search problem can be solved by a variety of quantum algorithms including a quantum exhaustive search in physical space or the coefficient space (Method IV), a quantum randomized search algorithm (Method V), or by employing a quantum hybrid algorithm (Method VI) depending on the nature of the problem. Quantum search algorithms offer a quadratic speed-up (in comparison to the traditional non-quantum approaches) and may well become methods of choice in many optimization applications once quantum computers become widely accessible and reach their true potential.