1 Introduction

1.1 Quantum computing

Quantum computing has provided novel approaches for solving computationally complex problems over the last decade by leveraging the inherent speedup(s) of quantum calculations compared to classical computing. Quantum superposition and entanglement are two key factors that give a massive speed up to calculations in the quantum domain compared to classical counterparts (Montanaro 2016; Jordan 2024; Horodecki et al. 2009). Because of this, addressing optimization problems by quantum computing is an appealing prospect. Multiple approaches, such as Grover’s algorithm (Grover 1996), adiabatic computation (AC) (Farhi et al. 2000), and quantum approximate optimization algorithm (QAOA) (Farhi et al. 2014), have been proposed to use quantum effects and, as such, have served as the basis for solving mathematically complex problems using quantum computing. The performance of classical algorithms has generally been found to be subpar when applied to larger dimensional problem spaces (National Academies of Sciences (2019)). On a multidimensional problem, classical machine learning optimization techniques frequently require a significant amount of CPU and GPU resources and long computation time. The reason for this is that ML techniques are needed to resolve NP-hard optimization problems (Dasari et al. 2020).

1.2 Vehicle routing problem

The vehicle routing problem is an intriguing optimization problem because of its many uses in routing and fleet management (Harwood et al. 2021), but its computational complexity is NP-hard (Kumar and Panneerselvam 2012; Mohanty et al. 2023). Moving automobiles as quickly and cheaply as feasible is always the objective. VRP has inspired a plethora of precise and heuristic approaches (Harwood et al. 2021; Srinivasan et al. 2018), all of which struggle to provide fast and trustworthy solutions. The VRP’s bare bones implementation comprises sending a single vehicle to deliver items to many client locations before returning to the depot to restock (Feld et al. 2019). By optimizing a collection of routes that are available and all initiate and conclude at a single node called the depot, the maximum reward sought by VRP is often expressed as the inverse of total distance traveled or mean service time. Even with just a few hundred delivery locations, finding the best solution to this problem is computationally challenging.

To be precise, in every VRP with parameters (n, k), there are (\(n-1\)) locations, k motor vehicles, and a depot D (Azad et al. 2023; Harwood et al. 2021). The solution is a collection of paths whereby each vehicle takes exactly one journey, and all k vehicles start and conclude at the same location, D. The best route is one that requires k vehicles to drive the fewest total miles. This problem may be thought of as a generalization of the well-known “traveling salesman” problem, whereby a set of k salesmen must service an aggregate of (\(n-1\)) sites with a single visit to each of those places being guaranteed (Harwood et al. 2021). In most practical settings, the VRP issue is complicated by other constraints, such as limited vehicle capacity or limited time for coverage. As a consequence, several other approaches, both classical and quantum, have been proposed as potential ways forwards. Currently, available quantum approaches for optimizing a system include the quantum approximate optimization algorithm (QAOA) (Azad et al. 2023), the quadratic unconstrained binary optimization (QUBO) (Glover et al. 2020), and quantum annealing (Irie et al. 2019; Crispin and Syrichas 2013; Office 2019).

1.3 Quantum support vector machine (QSVM)

The goal of the support vector machine (SVM) technique is to find the best line (or decision boundary) between two classes in n-dimensional space so that new data may be classified quickly. This optimum decision boundary is referred to as a hyperplane. The most extreme vectors and points that help construct the hyperplane are selected using SVM. The SVM method is based on support vectors, which are used to represent these extreme instances. Typically, a hyperplane cannot divide a data point in its original space. In order to find this hyperplane, a nonlinear transformation is applied to the data as a function. A feature map is a function that transforms the features of provided data into the inner product of data points, also known as the kernel (Havlicek et al. 2019; Rebentrost et al. 2014; Kariya and Behera 2024).

Quantum computing produces implicit calculations in high-dimensional Hilbert spaces using kernel techniques by physically manipulating quantum systems. Feature vectors for SVM in the quantum realm are represented by density operators, which are themselves encodings of quantum states. The kernel of a quantum support vector machine (QSVM) is made up of the fidelities between different feature vectors, as opposed to a classical SVM; the kernel conducts an encoding of classical input into quantum states (Havlicek et al. 2019; Leporini and Pastorello 2021).

1.4 Novelty and contribution

  • In this work, we propose a new method to solve the VRP using a machine learning approach through the use of QSVM.

  • In this context, we came across recent and older works in QSVM (Kariya and Behera 2024; Gentinetta et al. 2024; Rebentrost et al. 2014) and VQE algorithms (Cerezo et al. 2021), which are used to solve optimization problems such as VRP. However, none of them use a hybrid approach to arrive at a solution.

  • Our work implements this new approach of solving VRP using gate-model simulation of a 3-city or 4-city problem on a 6-qubit or 12-qubit system, respectively, using a parameterized circuit that is proposed as a solution to VRP.

  • We apply quantum encoding techniques such as amplitude encoding, angle encoding, higher-order encoding, IQP Encoding, and quantum algorithms such as QSVM, VQE, and QAOA to construct circuits for VRP and assess the results and summarize our findings.

  • We evaluate our solution using a variety of classical optimizers, as well as fixed and variable Hamiltonians to draw statistical conclusions.

1.5 Organization

The paper is organized as follows. Section 2 discusses the fundamental mathematical concepts such as QAOA, the Ising model, quantum support vector machine, amplitude encoding, angle encoding, higher-order encoding, IQP encoding, and VQE. Section 3 discusses the formulation and solution of VRP utilizing the concepts covered in the preceding section. Section 3.2 includes the fundamental components of circuits required to solve VRP using QSVM. Section 4 covers the outcomes of the QSVM simulation consisting of two sub-sections. Section 4.1 covers the outcome of simulation results of all the encoding schemes used, Finally, in Sect. 4.2, we conclude by comparing the results of QSVM solutions using various optimizers in the Qiskit platform on the VRP circuit and discuss the feasibility of higher qubit solutions as the future directions of research.

2 Background

Dealing with methods and processes for resolving combinatorial optimization problems is the foundation of solving routing challenges. The objective function is then derived by converting the mathematical models into their quantum counterparts. We arrive at the objective function’s solution by iteratively maximizing or minimizing the mathematical model. In this section, we provide an outline of our solution strategy’s key concepts.

2.1 QAOA

The quantum approximate optimization algorithm (QAOA) was proposed by Farhi et al. (2000, 2014) using an adiabatic quantum computation framework as the algorithm’s foundation. It is a hybrid algorithm because both classical and quantum approaches are utilized. Quantum adiabatic computation entails transitioning between the eigenstates of the driver Hamiltonian and the problem Hamiltonian. The Hamiltonian problem can be expressed as

$$\begin{aligned} C|z\rangle =\sum ^{m}_{\alpha =1}C_{\alpha }|z\rangle . \end{aligned}$$
(1)

It is well known that the combinatorial optimization problem may be effectively addressed by determining the eigenstate of C with the maximum energy. Likewise, we use driver Hamiltonian as

$$\begin{aligned} B=\sum ^n_{j=1}{\ }{\sigma }^x_j, \end{aligned}$$
(2)

Here, \({\sigma }^x_j\) denotes the \({\sigma }^x\) Pauli operator on bit \(z_j\), and B is the mixing operator. Let us additionally define two operators

$$\begin{aligned} U_C\left( \gamma \right) {=}e^{{-}i\gamma C}, U_B\left( \beta \right) {=}e^{{-}i\beta {B}}. \end{aligned}$$
(3)

This enables the system to develop under C for \(\gamma \) time and B for \(\beta \) time. Essentially, QAOA creates the following state,

$$\begin{aligned} |\varvec{\beta },\varvec{\gamma }\rangle = U_B\left( \beta _p\right) U_C\left( \gamma _p\right) \cdots U_B\left( \beta _2\right) U_C\left( \gamma _2\right) U_B\left( \beta _1\right) U_C\left( \gamma _1\right) |s\rangle \nonumber \\ \end{aligned}$$
(4)

Here, \(|s\rangle \) represents the superposition state of all inputs. The expectation value of the cost function \(\sum ^m_{\alpha {=1}}{\langle }\beta {,}\gamma ~{|C_\alpha }\)\({|}\beta {,}\gamma ~{\rangle }\) provides the solution, or an approximate solution to the problem (Zhou et al. 2020).

2.2 Ising model

The Ising model is a well-known mathematical representation of ferromagnetism in statistical mechanics (Singh 2020; Brush 1967). In the model, discrete variables (\(+1\) or \(-1\)) represent magnetic dipole moments of “spins” in one of two possible states. Each spin can interact with its neighbors because they are organized in a network, commonly a lattice (when there is periodic repetition in all directions of the local structure). The spins interact in pairs, with one value of energy when the two spins are identical and another value when they are dissimilar. However, heat reverses this tendency, permitting the formation of alternative structural phases. The model is a condensed representation of reality that allows phase transitions to be identified. The subsequent Hamiltonian describes the entire spin energy:

$$\begin{aligned} H_c=-\sum _{\left\langle i,j\right\rangle }{\ }J_{ij}{{\sigma }}_{{i}}{{\sigma }}_{{j}}-h\sum {{\sigma }}_{{i}}, \end{aligned}$$
(5)

where \(J_{ij}\) is the interaction between adjacent spins i and j, and h is an external magnetic field. At \(h=0\), the ground state is ferromagnetic if J is positive. At \(h=0\), the ground state is antiferromagnetic if J is negative in a bipartite lattice. Thus, for the purpose of clarity and within the scope of this paper, the Hamiltonian can be expressed as

$$\begin{aligned} H_c=-\sum _{\langle i,j\rangle }{\ }J_{ij}{\sigma }_{i}^z{\sigma }_{j}^z-\sum h_i{\sigma }_{i}^x. \end{aligned}$$
(6)

Here, \({\sigma }_z\) and \({\sigma }_x\) represent the z and x Pauli operators, respectively. For the sake of simplification, we can presume the following conditions to be ferromagnetic (\(J_{ij}>0\)) if there is no external impact on the spin: \(h=0\). The Hamiltonian may therefore be reformulated as follows:

$$\begin{aligned} H_c=-\sum _{\langle i,j\rangle }{\ }J_{ij}{\sigma }_{i}^z{\sigma }_{j}^z=-\sum _{\langle i,j\rangle }{\ }\sigma _{i}^z\sigma _{j}^z. \end{aligned}$$
(7)

2.3 Quantum support vector machine

SVM (Rebentrost et al. 2014; Kariya and Behera 2024) is a supervised algorithm that constructs hyperplane with \(\textbf{w} \cdot \textbf{x} + b = 0\) such that \(\textbf{w}\cdot \textbf{x} + b \ge 1\) for a training point \(\textbf{x}_i\) in the positive class, and \(\textbf{w}\cdot \textbf{x} + b \le -1 \) for a training point \(\textbf{x}_i\) in the negative class. During the training process, the algorithm aims to maximize the gap between the two classes, which is intuitive as we want to separate two classes as far as possible, in order to get a sharper estimate for the classification result of new data samples like \(\mathbf {x_0}\). Mathematically, we can see the objective of SVM is to find a hyperplane that maximizes the distance \(2/ |\textbf{w}|\) constraint to \(\mathbf {y_i}(\textbf{w}\cdot \mathbf {x_i} + b) \ge 1\). The normal vector \(\textbf{w}\) can be written as \(\textbf{w} = \sum _{i=1}^{M} \alpha _i \textbf{x}_i\) where \(\alpha _i\) is the weight of the \(i^{th}\) training vector \(\textbf{x}_i\). Thus, obtaining optimal parameters b and \(\alpha _i\) is the same as finding the optimal hyperplane. To classify the new vector is analogous to knowing which side of the hyperplane it lies, i.e., \( y_i(\textbf{x}_0) = sign(\textbf{w}.\textbf{x} + b) \). After having the optimal parameters, classification now becomes a linear operation. From the least-squares approximation of SVM, the optimal parameters can be obtained by solving a linear equation

$$\begin{aligned} \textbf{F}(b,\alpha _1,\alpha _2,\alpha _3,...,\alpha _M)^T = (0,y_1,y_2,y_3,...y_M)^T. \end{aligned}$$
(8)

In a general form of F, we adopt the linear kernels \(K_{i,j} = \kappa (\textbf{x}_i,\textbf{x}_j) = \textbf{x}_i. \textbf{x}_j \). Thus, to find the hyperplane parameters, we use matrix inversion of F: \((b,{\mathbf {\alpha }_i}^T)^T = {\tilde{F}^{ - 1}} (0,{\textbf{y}_i}^T)^T\).

2.3.1 Quantum kernels

The main inspiration of a quantum support vector machine is to consider quantum feature maps that lead to quantum kernel functions, which are hard to simulate in classical computers. In this case, the quantum computer is only used to estimate a quantum kernel function, which can be later used in kernel-based algorithms. For simplicity, assuming the data points \(x,z \in \mathcal {X}\), the nonlinear feature map of any data point is

$$\begin{aligned} \Phi (\varvec{x})=U(\varvec{x})\left| 0^n\right\rangle \left\langle 0^n\right| U^{\dagger }(\varvec{x}). \end{aligned}$$
(9)

The kernel function \(\kappa (x, z)\) can be computed as

$$\begin{aligned} \kappa (x, y)=|\langle \phi (x) \mid \phi (z)\rangle |^2. \end{aligned}$$
(10)

The state \(|\phi (x)\rangle \) can be prepared by using a unitary gate U(x), and thus \(|\phi (x)\rangle =U(x)|0\rangle \). Thus, the kernel fuction becomes

$$\begin{aligned} \kappa (x, z)=\left| \left\langle 0\left| U^{\dagger }(x) U(z)\right| 0\right\rangle \right| ^2. \end{aligned}$$
(11)

From the above, we can say that the kernel \(\kappa (x, z)\) is simply the probability of getting an all-zero string when the circuit \(U^{\dagger }(x) U(z)|0\rangle \) is measured or this kernel is an \(\left| 0^n\right\rangle \) to \(\left| 0^n\right\rangle \) transition probability of a particular unitary quantum circuit on n qubits (Havlicek et al. 2019; Glick et al. 2024). This can be implemented using the following kernel estimation circuit (Fig. 1).

Fig. 1
figure 1

Schematic diagram depicting quantum circuit for kernel estimation

Full size image

2.4 Amplitude encoding (AE)

In the process of amplitude embedding (Araujo et al. 2021), data is encoded into the amplitudes of a quantum state. A N-dimensional classical data point x is represented by the amplitudes of an n-qubit quantum state \(| \psi _x \rangle \) as

$$\begin{aligned} \left| \psi _x\right\rangle =\sum _{i=1}^N x_i|i\rangle \end{aligned}$$
(12)

where \(N=2^n\), \(x_i\) is the i-th element of x, and \(|i\rangle \) is the i-th computational basis state.

In order to encode any data point x into an amplitude-encoded state, we must normalize the same by following

$$\begin{aligned} | \psi _{x_{norm}} \rangle =\frac{1}{x_{norm}}\sum _{i=1}^N x_i|i\rangle , \end{aligned}$$
(13)

where \(x_{norm}\)= \(\sqrt{\sum _{i=1 }^N|x_i|^2}\).

2.5 Angle encoding (AgE)

While the above-described amplitude encoding expands into a complicated quantum circuit with huge depths, the angle encoding employs N qubits and a quantum circuit with fixed depth, making it favorable to NISQ computers (LaRose and Coyle 2020; Qiskit 2024). We define angle encoding as a method of classical information encoding that employs rotation gates (the rotation could be chosen along x, y, or z axis). In our scenario, the classical information consists of the node and edge weights assigned to the vehicle’s nodes and pathways which are further assigned as parameters to ansatz.

$$\begin{aligned} |\textbf{x}\rangle =\bigotimes _i^n R\left( \textbf{x}_i\right) \left| 0^n\right\rangle , \end{aligned}$$
(14)

where \({x}_i\) represents the classical information stored on the angle parameter of rotation operator R.

2.6 Higher-order encoding (HO)

Higher-order encoding is a variation of angle encoding where we have an entangled layer and an additional sequential operation of rotation angles of two entangled qubits (Qiskit 2024). This can be loosely defined as follows:

$$\begin{aligned} |\textbf{x}\rangle =\bigotimes _{i=2}^{n}R_(x_{i-1}.x_i)\bigotimes _{i=2}^{n-1}CX_{i,i+1}\bigotimes _{i=1}^n R\left( x_i\right) \left| 0^n\right\rangle . \end{aligned}$$
(15)

Similar to angle encoding, we are free to choose the rotation.

2.7 IQP encoding (IqpE)

IQP-style encoding is a relatively complicated encoding strategy. We encode classical information (Paddle Quantum 2024)

$$\begin{aligned} |{x}\rangle =\left( \textrm{U}_{\textrm{Z}}({x}) \textrm{H}^{\otimes n}\right) ^r\left| 0^n\right\rangle , \end{aligned}$$
(16)

where r is the depth of the circuit, indicating the repeating times of \(\textrm{U}_{\textrm{Z}}({x}) \textrm{H}^{\otimes n}\). \(\textrm{H}^{\otimes n}\) is a layer of Hadamard gates acting on all qubits. \(\textrm{U}_{\textrm{Z}}(\textbf{x})\) is the key step in IQP encoding scheme:

$$\begin{aligned} \textrm{U}_{\textrm{Z}}(\textrm{x})=\prod _{[i, j] \in S} R_{Z_t Z_j}\left( x_i x_j\right) \bigotimes _{k=1}^n R_z\left( x_k\right) , \end{aligned}$$
(17)

where S is the set containing all pairs of qubits to be entangled using \(R_{Z Z}\) gates. First, we consider a simple two-qubit gate: \(R_{Z_1 Z_2}(\theta )\). Its mathematical form \(e^{-i \frac{i}{2} Z_1 \otimes Z_2}\) can be seen as a two-qubit rotation gate around ZZ, which makes these two qubits entangled.

2.8 VQE

Variational quantum eigensolver (VQE) is another hybrid quantum-classical algorithm used for the estimation of the eigenvalue of a matrix or Hamiltonian H (Peruzzo et al. 2014) of significant size. The primary objective of this approach is to ascertain a trial qubit state from a wave function \(| \psi (\mathbf {\theta }) \rangle \) that relies on a collection of parameters \(\mathbf {\theta }= \theta _1,\theta _2,\cdots \), which are often referred to as the variational parameters. The expectation of an observable or Hamiltonian H in a state \(| \psi (\mathbf {\theta }) \rangle \) can be expressed as follows:

$$\begin{aligned} E(\mathbf {\theta })= \langle \psi (\mathbf {\theta }) | H | \psi (\mathbf {\theta }) \rangle . \end{aligned}$$
(18)

By spectral decomposition

$$\begin{aligned} H=\lambda _{1}| \psi \rangle _{1}\langle \psi |_{1}+\lambda _{2}| \psi \rangle _{2}\langle \psi |_{2}+\ldots +\lambda _{n}| \psi \rangle _{n}\langle \psi |_{n},\nonumber \\ \end{aligned}$$
(19)

where \(\lambda _i\) and \({| \psi \rangle }_i\) are the matrix H’s eigenvalues and eigenstates, respectively. Additionally, because H’s eigenstates are orthogonal, \(\left\langle \psi _{i} \mid \psi _{j}\right\rangle =0\) If \(i \ne j\). The wave function \(| \psi (\mathbf {\theta }) \rangle \) can be expressed as a superposition of eigenstates.

$$\begin{aligned} | \psi (\mathbf {\theta }) \rangle =\alpha _{1}(\mathbf {\theta })| \psi \rangle _{1}+\alpha _{2}(\mathbf {\theta })| \psi \rangle _{2}+\ldots +\alpha _{n}(\mathbf {\theta })| \psi \rangle _{n}. \end{aligned}$$
(20)

Hence, the expectation is given by

$$\begin{aligned} E(\mathbf {\theta })= & {} |\alpha _{1}(\mathbf {\theta }) |^2 \lambda _{1}+|\alpha _{2} (\mathbf {\theta }) |^2\lambda _2+\ldots +|\alpha _{n}(\mathbf {\theta }) |^2 \lambda _{n}.\nonumber \\ \end{aligned}$$
(21)

Hence, \( E(\mathbf {\theta }) \ge \lambda _{\min }\). The VQE method involves the iterative adjustment of the parameters \(\mathbf {\theta }=\theta _{1}, \theta _{2}, \ldots \) in order to minimize the value of \(E(\mathbf {\theta })\). This property of VQE is advantageous when attempting to solve combinatorial optimization problems; specifically, the approach involves using a parameterized circuit to establish the trial state of the algorithm, with the cost function denoted as \(E(\mathbf {\theta })\), which is also the expected value of the Hamiltonian in this state. It is possible to derive the ground state of the desired Hamiltonian by iteratively minimizing the cost function. A classical optimizer uses a quantum computer to calculate its gradient and assess the cost function at each step of the optimization process.

3 Methodology

3.1 Modelling VRP in QSVM

By mapping the cost function to an Ising Hamiltonian \(H_c\), the vehicle routing problem can be solved (Lucas 2014). The solution to the problem is determined by minimizing the Ising Hamiltonian \(H_c\). Consider a graph with n vertices and \(n-1\) edges and an arbitrary connectivity. Assuming we must route a vehicle between two non-adjacent vertices in the graph, consider a binary decision variable \(x_{ij}\) whose value is 1 if there is an edge between i and j with an edge weight \(w_{ij}>0\) and 0 otherwise. Now, the VRP problem necessitates \(n \times (n-1)\) selection variables. We define two sets of nodes for each edge from \(i\rightarrow j\): source \(s\left[ i\right] \) and target t[j]. \(s\left[ i\right] \) contains the nodes j to which i sends an edge \(j\ \epsilon \ s[i]\). The collection \(t\left[ j\right] \) comprises the nodes i to which the node i delivers the edge \(i\ \epsilon \ t[j]\). The VRP is defined as follows (Azad et al. 2023, Qiskit: Vehicle Routing 2024):

$$\begin{aligned} VRP(n,k)=\mathop {min}_{{\left\{ x_{ij}\right\} }_{i\rightarrow j}\in \{0,1\}}\ \sum _{i\rightarrow j}{\ }w_{ij}x_{ij}, \end{aligned}$$
(22)

where k and n represent the number of vehicles and locations respectively; there are \(n-1\) locations for vehicles to traverse if the starting point is considered to be the 0th location or Depot D. Notably, the following restrictions apply to this (Mohanty et al. 2023):

$$\begin{aligned} \sum _{j\in ~{s}~[i]}{\ }x_{ij}= & {} 1,{\forall }i\in \{1,\cdots ,n-1\}, \nonumber \\ \sum _{j\in ~{t}~[i]}{\ }x_{ji}= & {} 1,{\forall }i\in \{1,\cdots ,n-1\}, \nonumber \\ \sum _{{j}{\in }~{s}~{[}{0}{]}}{{\ }}{x}_{0j}= & {} k, \nonumber \\ \sum _{j\in ~{t}~[0]}{\ }x_{j0}= & {} k \nonumber \\ u_i-u_j+Q x_{i j}\le & {} Q-q_j, \forall i \sim j, i, j \ne 0, \nonumber \\ q_i \le u_i\le & {} Q, \forall i, i \ne 0. \end{aligned}$$
(23)

The first two restrictions establish that each node may only be visited once by the delivering vehicle. The middle two limitations enforce the requirement that after product delivery, the vehicle must return to the depot. The last two restrictions enforce the requirements for eliminating sub-tours and are constrained on \(u_i\), with \(Q>q_j>0\), and \(u_i,Q, q_i \in \mathbb {R}\). For the VRP equation and constraints, the VRP Hamiltonian can be expressed as follows (Azad et al. 2023).

$$\begin{aligned} H_{VRP}= & {} H_a+H_b+H_c+H_d+H_e, \nonumber \\ H_a= & {} ~\sum _{i~\rightarrow j}{w_{ij}x_{ij}},\nonumber \\ H_b= & {} A\sum _{i\in 1,\cdots ,n-1}{\ }{\left( 1-\sum _{j\in ~{s}~[i]}{\ }x_{ij}\right) }^2, \nonumber \\ H_c= & {} A\sum _{i\in 1,\cdots ,n-1}{\ }{\left( 1-\sum _{j\in ~{t}[i]}{\ }x_{ji}\right) }^2, \nonumber \\ H_d= & {} A{\left( k-\sum _{j\in ~{s}[0]}{\ }x_{0j}\right) }^2, \nonumber \\ H_e= & {} A{\left( k-\sum _{j\in ~{t}[0]}{\ }x_{j0}\right) }^2. \end{aligned}$$
(24)

\(A > 0\) is indicative of a constant. The vector representation of the collection of all binary decision variables \(x_{ij}\) is

$$\begin{aligned} \overrightarrow{\varvec{{x}}}={\left[ x_{(0,1)},x_{(0,2)},\cdots x_{(1,0)},x_{(1,2)},\cdots x_{(n-1,n-2)}\right] }^{\varvec{{T}}}. \end{aligned}$$
(25)

Using the antecedent vector, we can construct two new vectors for each node: \(\overrightarrow{z}_{S[i]}\) and \(\overrightarrow{z}_{T[i]}\) (at the start of the section, we defined two sets for source and target nodes, so two vectors will represent them).

$$\begin{aligned} \overrightarrow{z}_{S\left[ i\right] }= & {} \textbf{x} \ni x_{ij}=1,\ x_{kj}\!=\!0\,\ k\ne i\,\ \ \forall j,k\!\in \! \{0,\cdots ,n-1\}, \nonumber \\ \overrightarrow{z}_{T\left[ i\right] }= & {} \textbf{x} \ni x_{ji}=1,\ x_{jk}\!=\!0\,\ k\ne i\,\ \ \forall j,k\! \in \! \{0,\cdots ,n-1\}. \nonumber \\ \end{aligned}$$
(26)
$$\begin{aligned} \sum _{j \in \text{ s }[i]} x_{i j}= & {} \textbf{z}_{S[i]}^{\textrm{T}} \textbf{x}, \nonumber \\ \sum _{j \in \text{ t }[i]} x_{j i}= & {} \textbf{z}_{T[i]}^{\textrm{T}} \textbf{x}. \end{aligned}$$
(27)

These vectors will contribute to the development of the QUBO model of VRP (Date et al. 2019; Glover et al. 2020; Kochenberger et al. 2014; Guerreschi 2021). The QUBO model of a connected graph \(G=(N,V)\) is specified as follows:

$$\begin{aligned} f(x)_{QUBO}=\mathop {min}_{x\in \{0,1\}(N\times V)}x^TQx+g^Tx+c, \end{aligned}$$
(28)

where Q is a quadratic edge weight coefficient, g is a linear node weight coefficient, and c is a constant. To determine the coefficients in the QUBO formulations of \(H_{VRP}\) as shown in Eq. 24, the equations in Eq. 27 are first substituted in terms of \(H_{b}\) and \(H_{c}\), respectively. Subsequently, the expression of \(H_{VRP}\) is expanded and rearranged in accordance with Eq. 28.

$$\begin{aligned} H= & {} A \sum _{i=0}^{n-1}\left[ z_{S[i]} z_{S[i]}^T+z_{T[i]} z_{T[i]}^T\right] \textbf{x}^2 \nonumber \\{} & {} +w^T \textbf{x}-2 A \sum _{i=1}^{n-1}\left[ z_{S[i]}^T+z_{T[i]}^T\right] \textbf{x} \nonumber \\{} & {} - 2 A k\left[ z_{S[0]}^T+z_{T[0]}^T\right] \textbf{x}\!+2 A(n\!-1)\!+2 A k^2. \end{aligned}$$
(29)

Hence, in the QUBO formulation of the Eq. 24, the coefficients \(\textrm{Q}(n(n-1) \times n(n-1))\), \(\textrm{g}(n(n-1) \times 1)\), and \(\textrm{c}\) are derived. The coefficients associated with the QUBO formulation of Eq. 24 are shown below.

$$\begin{aligned} Q= & {} A\left[ \left[ z_{T[0]}, \ldots , z_{T[n-1]}\right] ^{T}\left[ z_{T[0]}, \ldots , z_{T[n-1]}\right] \right. \nonumber \\{} & {} \left. +\left( \mathbb {I}_{n} \otimes \mathbb {J}(n-1, n-1)\right) \right] , \nonumber \\ g= & {} W-2 A k\left( \left( e_{0} \otimes \mathbb {J}_{n-1}\right) +\left[ z_{T[0]}\right] ^{T}\right) , \nonumber \\ \quad{} & {} +2 A\left( \mathbb {J}_{n} \otimes \mathbb {J}_{n-1}\right) , \nonumber \\ c= & {} 2 A(n-1)+2 A k^{2}. \end{aligned}$$
(30)

\(\mathbb {J}\) is the matrix containing all ones, \(\mathbb {I}\) and \(e_0={\left[ 1,0,\cdots ..,0\right] }^T\) are the identity matrices. The binary decision variable \(x_{ij}\) is converted to the spin variable \(s_{ij} \in \left\{ -1,1\right\} \) using the formula \(x_{ij}=(s_{ij}+1)/2\).

Fig. 2
figure 2

a Circuit example illustrating gate operations for \({H}_{\textrm{c}}\). b Circuit example displaying gate selections with an additional u gate for \({H}_{\textrm{mxr}}\)

Full size image

From the aforementioned equations, we may expand Eq. 28 to form the Ising Hamiltonian of VRP (Glover et al. 2020).

$$\begin{aligned} H_{Is}=-\sum _i{\ }\sum _{i<j}{\ }J_{ij}s_is_j-\sum _i{\ }h_is_i+d. \end{aligned}$$
(31)

The following are definitions for the terms \(J_{ij}, h_i\), and d:

$$\begin{aligned} J_{ij}= & {} \ -\frac{Q_{ij}}{2},\ \forall \ i<j,\nonumber \\ h_i= & {} \frac{g_i}{2}+\sum {\frac{Q_{ij}}{4}+\ \sum {\frac{Q_{ji}}{4}\ }\ }, \nonumber \\ d= & {} c+\sum _i{\ }\frac{g_i}{2}+\sum _i{\ }\sum _j{\ }\frac{Q_{ij}}{4}. \end{aligned}$$
(32)

3.2 Analysis and circuit building

3.2.1 VRP

In the current section, we proceed to create a circuit based on gates using the IBM gate model. The implementation of this model is carried out using the Qiskit framework (et al. 2024), enabling us to effectively execute the aforementioned formulation. In the context of a given vehicle routing problem (VRP) that incorporates qubits, the initial state is established as \(| + \rangle ^{\otimes n(n-1)}\). This state represents the ground state of \(H_{mxr}\), which is achieved by applying the Hadamard gate to each qubit that has been initialized to the zero state. Subsequently, we proceed to build the subsequent state.

$$\begin{aligned} | \beta ,\gamma \rangle= & {} e^{-iH_{mxr}\beta _{p}}e^{-iH_{c}\gamma _{p}}...\nonumber \\{} & {} ...e^{-iH_{mxr}\beta _{0}}e^{-iH_{c}\gamma _{0}}| + \rangle ^{n\otimes (n-1)}. \end{aligned}$$
(33)

The energy E of the state \(| \beta ,\gamma \rangle \) is computed using the expectation of \(H_{c}\) from Eq. 18. Again, based on the Ising model, the term \(H_{c}\) can be expressed in terms of Pauli operators as follows:

$$\begin{aligned} {H}_{\textrm{c}}=-\sum _{i} \sum _{i<j} J_{ij}\sigma _i^z\sigma _j^z-\sum _{i} h_i\sigma _i^z-d. \end{aligned}$$
(34)

Thus, the expression for a single term of state in \(| \beta ,\gamma \rangle \) as \(\beta _0,\gamma _0\) reads as follows:

\(e^{-i{H}_{mxr}\beta _0}e^{-i{H}_{c}\gamma _0}.\) The first term \({H}_{\mathrm {c\ }}\) can be expanded as follows:

$$\begin{aligned} {e}^{iJ_{ij}\gamma _0\sigma _i\sigma _j}= & {} \cos J_{ij}\gamma _0I+i\ \sin J_{ij}\gamma _0\sigma _i\sigma _j,\nonumber \\= & {} {e}^{i{J}_{ij}{\gamma }_{0}}| 00 \rangle \langle 00 |+{e}^{-i{J}_{ij}{\gamma }_{0}}| 01 \rangle \langle 01 |\nonumber \\{} & {} +{e}^{-i{J}_{ij}{\gamma }_{0}} | 10 \rangle \langle 10 |+ {e}^{i{J}_{ij}{\gamma }_{0}}| 11 \rangle \langle 11 |, \nonumber \\= & {} M \end{aligned}$$
(35)

It can be noted that, by the application of CNOT(CX) gate before and after “M,” the diagonal elements of the above matrix can be swapped.

$$\begin{aligned} CX(M)CX= & {} {e}^{i{J}_{ij}{\gamma }_{0}}| 00 \rangle \langle 00 | - {e}^{-i{J}_{ij}{\gamma }_{0}}| 01 \rangle \langle 01 |\nonumber \\{} & {} +{e}^{-i{J}_{ij}{\gamma }_{0}} | 10 \rangle \langle 10 | - {e}^{i{J}_{ij}{\gamma }_{0}}| 11 \rangle \langle 11 |, \end{aligned}$$
(36)

Expanding the matrix and observing the upper and lower blocks, we can rewrite

$$\begin{aligned} \left[ \begin{matrix}1&{}0\\ 0&{}1\end{matrix}\right] \otimes \left[ \begin{matrix}{e}^{i{J}_{ij}{\gamma }_{0}}&{}0\\ 0&{}{e}^{-i{J}_{ij}{\gamma }_{0}}\end{matrix}\right] = I\otimes {e}^{i{J}_{ij}{\gamma }_{0}}\left[ \begin{matrix}1&{}0\\ 0&{}{e}^{-2i{J}_{ij}{\gamma }_{0}}\end{matrix}\right] .\nonumber \\ \end{aligned}$$
(37)

\(\left[ \begin{matrix}1&{}0\\ 0&{}{e}^{-2i{J}_{ij}{\gamma }_{0}}\end{matrix}\right] \) is a phase gate. Looking at the 2nd term of \({H}_{\textrm{c}}\), we get

$$\begin{aligned} {H}_{\textrm{c}}= & {} \sum _{i} h_i\sigma _{i}^{z}, \nonumber \\ e^{ih_{i}\sigma _i}= & {} {\cos h}_i{\gamma _oI+i\sin \gamma _0\sigma _i},\nonumber \\= & {} \cos h_i{\gamma _o\left[ \begin{matrix}1&{}0\\ 0&{}1\\ \end{matrix}\right] +i\sin h_i\gamma _0\left[ \begin{matrix}1&{}0\\ 0&{}-1\\ \end{matrix}\right] },\nonumber \\= & {} \left[ \begin{matrix}e^{ih_i\gamma _0}\ {} &{}0\\ 0&{}e^{-ih_i\gamma _0}\\ \end{matrix}\right] . \end{aligned}$$
(38)
Fig. 3
figure 3

Plot illustrating different encoding methods for two qubits. a Amplitude encoding, b angle encoding, c higher-order encoding, and d IQP encoding

Full size image

Figure 2a picturizes the fundamental circuit with two qubits and gate selections for \({H}_{\textrm{c}}\). Similarly, \(H_{mxr}\) is merely a rotation along the X axis, as depicted by the U gate in Fig. 2b.

The above sample circuits can be used for the solution of VRP combined with the VQE and QAOA approach. However, in this paper, we are focusing on a machine learning solution of VRP by use of QSVM; thus, we need to construct a QSVM circuit using various encoding schemes. Simple interpretation and implementation of encoding schemes are described in upcoming sub-sections.

3.2.2 Amplitude encoding

As we look into AE, a single qubit state is represented by

$$\begin{aligned} | \psi \rangle (\theta )=\cos (\theta /2)| 0 \rangle +\sin (\theta /2)| 1 \rangle , \end{aligned}$$
(39)

for two qubits

$$\begin{aligned} |\psi (\theta )\rangle= & {} \alpha |00\rangle +\beta |01\rangle +\gamma |10\rangle +\delta |11\rangle , \nonumber \\= & {} |0\rangle (\alpha |0\rangle +\beta |1\rangle )+|1\rangle (\gamma |0\rangle +\delta |1\rangle ), \nonumber \\= & {} |0\rangle \sqrt{\left( \alpha ^2+\beta ^2\right) }\left( \frac{\alpha |0\rangle +\beta |1\rangle }{\sqrt{\alpha ^2+\beta ^2}}\right) \nonumber \\{} & {} +|1\rangle \sqrt{\gamma ^2+\delta ^2} \frac{\gamma |0\rangle +\delta |1\rangle }{\sqrt{\gamma ^2+\delta ^2}}. \end{aligned}$$
(40)

Now, applying Ctrl U and Anti-CTRL U on the above state, we achieve

$$\begin{aligned}{} & {} |0\rangle \sqrt{\alpha ^2+\beta ^2}|0\rangle +|1\rangle \sqrt{\gamma ^2+\delta ^2}|0\rangle \nonumber \\= & {} \left( \sqrt{\alpha ^2+\beta ^2}|0\rangle +\sqrt{\gamma ^2+\delta ^2}|1\rangle \right) |0\rangle . \end{aligned}$$
(41)

Here, \(\theta _1 = \tan ^{-1} \frac{\sqrt{\gamma ^2+ \delta ^2}}{\sqrt{\alpha ^2 + \beta ^2}}\), \(\theta _2= \tan ^{-1}\frac{\delta }{\gamma }\), \(\theta _3= \tan ^{-1}\frac{\beta }{\alpha }\) Combining VRP and amplitude encoding circuit eliminates the need of Hadamard gates and \(H_{mxr}\) components, and we end up with the following skeleton circuits (Fig. 3a).

3.2.3 Angle encoding

For a two-qubit scenario, angle encoding translates to the following example. We define the \(R_y\) gate as follows:

$$\begin{aligned} R_y(\theta )= & {} e^{-i Y \theta / 2}=\cos \frac{\theta }{2}-i \sin {\theta / 2} Y, \nonumber \\= & {} \left[ \begin{array}{ll} \cos \theta / 2 &{} -\sin \theta / 2 \\ \sin \theta / 2 &{} \cos \theta / 2 \end{array}\right] . \end{aligned}$$
(42)
$$\begin{aligned}{} & {} |00\rangle \nonumber \\{} & {} \xrightarrow [R_y(\theta _2)]{R_y(\theta _1)}\left( \cos \frac{\theta _1}{2}|0\rangle +\sin \frac{\theta _1}{2}|1\rangle \right) \left( \cos \frac{\theta _2}{2}|0\rangle +\sin \frac{\theta _2}{2}|1\rangle \right) , \nonumber \\= & {} \cos \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|00\rangle +\cos \frac{\theta _1}{2} \cdot \sin \frac{\theta _2}{2}|01\rangle \nonumber \\{} & {} +\sin \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|10\rangle +\sin \frac{\theta _1}{2} \sin \frac{\theta _2}{2}|11\rangle \nonumber \\{} & {} \xrightarrow []{CNOT} \cos \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|00\rangle +\cos \frac{\theta _1}{2} \cdot \sin \frac{\theta _2}{2}|01\rangle \nonumber \\{} & {} +\sin \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|11\rangle +\sin \frac{\theta _1}{2} \sin \frac{\theta _2}{2}|10\rangle . \end{aligned}$$
(43)

3.2.4 Higher-order encoding

For a two-qubit scenario, HO encoding translates to the following: We define the \(R_y\) gate as follows:

$$\begin{aligned} R_y(\theta )= & {} e^{-i Y \theta / 2}=\cos \frac{\theta }{2}-i \sin {\theta / 2} Y, \nonumber \\= & {} \left[ \begin{array}{ll} \cos \theta / 2 &{} -\sin \theta / 2 \\ \sin \theta / 2 &{} \cos \theta / 2 \end{array}\right] . \end{aligned}$$
(44)
$$\begin{aligned}{} & {} |00\rangle \nonumber \\{} & {} \xrightarrow [R_y(\theta _2)]{R_y(\theta _1)}\left( \cos \frac{\theta _1}{2}|0\rangle +\sin \frac{\theta _1}{2}|1\rangle \right) \left( \cos \frac{\theta _2}{2}|0\rangle +\sin \frac{\theta _2}{2}|1\rangle \right) , \nonumber \\= & {} \cos \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|00\rangle +\cos \frac{\theta _1}{2} \cdot \sin \frac{\theta _2}{2}|01\rangle \nonumber \\{} & {} +\sin \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|10\rangle +\sin \frac{\theta _1}{2} \sin \frac{\theta _2}{2}|11\rangle \nonumber \\{} & {} \xrightarrow [R_y(\theta _1.\theta _2)]{CNOT} \cos \frac{\theta _1}{2} \cdot \cos \frac{\theta _2}{2}|0\rangle \left( \cos \frac{ \theta _1 \cdot \theta _2}{2}|0\rangle +\sin \frac{\theta _1 \cdot \theta _2}{2}|1\rangle \right) \nonumber \\{} & {} +\cos \frac{\theta _1}{2} \cdot \sin \frac{\theta _2}{2}|0\rangle \left( -\sin \frac{ \theta _1 \cdot \theta _2}{2}|0\rangle +\cos \frac{\theta _1 \cdot \theta _2}{2}|1\rangle \right) \nonumber \\{} & {} +\sin \frac{\theta _1}{2} \cos \frac{\theta _2}{2}|1\rangle \left( \cos \frac{\theta _1 \cdot \theta _2}{2}|0\rangle +\sin \frac{\theta _1 \cdot \theta _2}{2}|1\rangle \right) \nonumber \\{} & {} +\sin \frac{\theta _1}{2} \cdot \sin \frac{\theta _2}{2}|1\rangle \left( -\sin \frac{\theta _1 \cdot \theta _2}{2} |0\rangle +\cos \frac{\theta _1 \cdot \theta _2}{2}|1\rangle \right) . \nonumber \\ \end{aligned}$$
(45)

3.2.5 IQP encoding

For a two-qubit scenario, IqpE translates to the following:

$$\begin{aligned}{} & {} |00\rangle \xrightarrow {H_1 H_2} |++\rangle , \nonumber \\= & {} \frac{1}{2}(|00\rangle +|01\rangle +|10\rangle +|11\rangle ), \nonumber \\{} & {} \xrightarrow [R_Z(\theta _2)]{R_Z(\theta _1)} \frac{1}{2}\left( |00\rangle +e^{i \theta _2}|01\rangle +e^{i \theta _1}|10\rangle +e^{i\left( \theta _1+\theta _2\right) }|11\rangle \right) \nonumber \\{} & {} \left. \overset{ \text{ CNOT } }{\longrightarrow } \frac{1}{2}(|00\rangle +e^{i \theta _2}|01\rangle +e^{i \theta _1}\left| 11\right\rangle +e^{i\left( \theta _1+\theta _2\right) }|10\rangle \right) \nonumber \\{} & {} \xrightarrow {R_Z(\theta _1.\theta _2)} \frac{1}{2}(|00\rangle +e^{i \theta _2(1+\theta _1)}|01\rangle +e^{i \theta _1(1+\theta _2)}|11\rangle \nonumber \\{} & {} +e^{i \theta _1(1+\theta _2)}|10\rangle +e^{i(\theta _1+\theta _2)}|11\rangle ). \end{aligned}$$
(46)

4 Results

4.1 VQE simulation of QSVM and VRP

We build the Hamiltonian with a uniform distribution of weights between 0 and 1, and then run it along with the ansatz via IBM’s three available VQE optimizers (COBYLA, L_BFGS_B, and SLSQP). We run the circuit up to two layers and gather data using all of the available optimizers. We run the experiment again with a fixed Hamiltonian and, subsequently, a set of variable Hamiltonians to see whether the QSVM and encoding approach can effectively reach the classical minimum. Our results indicate that COBYLA is the most efficient optimizer, followed by SLSQP and L BFGS B. In the sections that follow, we will have a look at the results obtained using various QSVM encoding schemes. We define two terms—accuracy and error—in the context of outcomes’ interpretability. An error occurs when the solution deviates from the classical minimum more often than it reaches it, whereas accuracy is defined as the number of times the solution reaches the classical minimum. Percentages based on the distribution of the outcomes are used to evaluate both terms.

$$\begin{aligned} Acc= & {} \frac{N}{T}, \nonumber \\ Err= & {} \frac{T - N}{T}. \end{aligned}$$
(47)

\(T =\) Total number of Simulation runs

\(N =\) Total number of times solution reaches classical minimum

4.1.1 Amplitude encoding

With a large number of gates, the AE circuit has proven to be the most complex of all encoding circuits. We can simulate no more than six qubit computations due to this complexity. Despite its complexity, AE has a high, nearly perfect accuracy rate (\(100\%\)) and a very low error rate (\(0\%\)) for 50-iteration fixed Hamiltonian simulations. The trend is present in both the first and second layers. The first layer accuracy for a variable Hamiltonian simulation is \(96\%\), and the second layer accuracy is \(94\%\) across all optimizers. Figure 4 depicts the results of 50 iterations of simulating SVM with amplitude encoding on a VRP circuit with fixed and variable Hamiltonian. The decline in accuracy, however, can be attributed to simulation or computational errors, as all the errors are greater than 100% and are therefore considered aberrations. Most likely, the simulation hardware cannot accommodate the VQE procedure.

Fig. 4
figure 4

Plot illustrating Amplitude encoding results for QSVM solution of VRP. a Amplitude encoding 6 qubits fix Hamiltonian and b amplitude encoding 6 qubits variable hamiltonian

Full size image
Fig. 5
figure 5

Plot illustrating angle encoding results for QSVM solution of VRP. a Angle encoding 6 qubits fix Hamiltonian, b angle encoding 12 qubits fix Hamiltonian, c angle encoding 6 qubits variable Hamiltonian, and d angle encoding 12 qubits variable Hamiltonian

Full size image
Fig. 6
figure 6

Plot illustrating Higherorder encoding results for QSVM solution of VRP. a Higher-order encoding 6 qubits fix Hamiltonian, b higher-order encoding 12 qubits fix Hamiltonian, c higher-order encoding 6 qubits variable Hamiltonian, d higher-order encoding 12 qubits variable Hamiltonian

Full size image
Fig. 7
figure 7

Plot illustrating IQP encoding results for QSVM solution of VRP. a IQP encoding 6 qubits fix Hamiltonian, b IQP encoding 12 qubits fix Hamiltonian, c IQP encoding 6 qubits variable Hamiltonian, d IQP encoding 12 qubits variable Hamiltonian

Full size image

4.1.2 Angle encoding

Angle encoding is the second encoding (Fig. 5), following amplitude encoding; we have experimented with SVM VRP simulation, which yields high accuracy and low error rates. Observing Tables 1 and 2, angle encoding is the second most precise encoding employed in our investigations. For fixed Hamiltonian simulations over 50 iterations with 6 qubits angle encoding, the first layer, including all optimizers, achieves \(100\%\) accuracy and \(0\%\) error. In the second layer simulation (over 50 iterations), the accuracy decreases to \(98\%\) for COBYLA, \(96\%\) for SLSQP, and \(86\%\) for L_BGFS_B, which is a greater decrease than the other two. These declines are attributable to optimizer-dependent statistical errors. Similarly, for 12-qubit simulations of SVM VRP, the accuracy rates are higher in the first layer, which consists of COBYLA at \(100\%\), SLSQP at \(92\%\), and L_BGFS_B at \(88\%\), reiterating that the accuracy is highly dependent on the optimizer. As we move to the second layer of 12-qubit simulations on fixed Hamiltonian, we observe a decline in precision as the level of optimization rises. In this case, COBYLA winds up with \(80\%\), L_BGFS_B with \(70\%\), and SLSQP with \(84\%\). Here, SLSQP’s accuracy loss is less than that of the other two optimizers. The variable Hamiltonian with 12 qubits demonstrates a comparable trend. On the initial layer, we observe high accuracy with COBYLA at \(96\%\), L_BGFS_B at \(86\%\), and SLSQP at \(90\%\). Moving to the second stratum, the accuracy figures drop significantly, with COBYLA at \(76\%\) and L_BGFS_B at \(62\%\), while SLSQP maintains excellent accuracy at \(86\%\). In every scenario of our investigation, it is evident that over optimization reduces accuracy rates.

4.1.3 Higher-order encoding

After amplitude and angle encoding, higher-order encoding is the third most prevalent encoding in our SVM VRP simulation experiment. This is also the third most accurate encoding in our experiment. For both 6 qubit and 12 qubit simulations, HO encoding yields moderately accurate results; however, as the number of circuit layers is increased, the accuracy of the HO encoding scheme deteriorates, rendering it inappropriate. Figure 6 depicts the statistics of the HO encoding scheme for fixed and variable Hamiltonian simulations of SVM VRP circuits over 50 iterations for both 6 qubit and 12 qubit simulations. COBYLA achieves \(78\%\) accuracy for a 6-qubit HO encoding circuit on a fixed Hamiltonian, while L_BGFS_B achieves \(66\%\) accuracy and SLSQP achieves \(70\%\) accuracy. As we proceed to the second layer, the accuracy considerably decreases, with COBYLA at \(34\%\) and SLSQP, L_BGFS_B at \(16\%\), respectively. Similar trends can be observed in variable Hamiltonian simulations of HO encoding with 6 qubits, with COBYLA at \(76\%\), SLSQP at \(62\%\), and L_BGFS_B at \(58\%\) for the first layer; for the second layer, the accuracy drops to \(36\%\), \(34\%\), and \(36\%\) for COBYLA, L_BGFS_B, and SLSQP, respectively. The 12-qubit simulation yields superior results than the 6-qubit simulation and improves COBYLA’s accuracy. For fixed Hamiltonian simulations, COBYLA achieves an accuracy of \(92\%\), compared to \(78\%\) for 6 qubit. For variable Hamiltonian simulations, COBYLA stores \(76\%\) for 6 qubit in the first layer, and \(92\%\) for 12 qubit in the first layer. The tendencies for L_BGFS_B and SLSQP are ambiguous for both cases (fixed and variable Hamiltonian simulations); it is reassuring to conclude that an increase in layer decreases accuracy and that COBYLA outperforms the other two optimizers and ensures stable performance.

4.1.4 IQP encoding

IQP encoding (Fig. 7) is the last and least accurate encoding in our experiment to simulate an SVM VRP circuit. The results are plotted in Fig. 7 and in Tables 1 and 2. As we can see from the figures and tables, accuracy is consistently poor for fixed and variable Hamiltonian simulations in both 6-qubit and 12-qubit circuits. The accuracy further declines as layers increase. Hence, this encoding is unsuitable in our experiment of SVM VRP circuits.

Table 1 For 6- and 12-qubit VRP circuits using SVM with 2 layers, this shows the accuracy and error with reference to classical minimum (over 50 iterations) for VQE simulations over a fixed Hamiltonian; utilizing amplitude, angle, higher-order, and IQP encoding schemes, over the use of COBYLA, SLSQP and L_BGFS_B optimizers
Full size table
Table 2 For 6- and 12-qubit VRP circuits using SVM on with two layers, this shows the accuracy and error with reference to classical minimum (over 50 iterations) for VQE simulations on variable Hamiltonians utilizing amplitude, angle, higher-order, and IQP encoding schemes, over the use of COBYLA, SLSQP and L_BGFS_B optimizers
Full size table

4.2 Inferences from simulation

As we scan through the results of SVM VRP simulations across the encoding schemes, we observe some clear and distinct trends regarding the experiment. Tables 1 and 2 summarize the results obtained from the plots of all the encoding schemes used in this experiment. We list these trends as our outcomes of this experiment in the below points:

  • The approach to solving VRP using machine learning is successful and is capable of accomplishing the same or a superior result than the conventional approach using VQE and QAOA.

  • The use of encoding/decoding schemes can serve the purpose of creating superposition and entanglement and eliminate the additional effort required to construct the mixer Hamiltonian when solving the VRP using the standard approach of QAOA and VQE.

  • While the standard approach to solving VRP or any combinatorial optimization problem requires a few layers of circuit depth (2 in most cases), we are able to achieve the same on the first layer itself with this approach, proving that it is more efficient than the standard approach.

  • We also observe a distinct trend that as the number of layers increases, the accuracy decreases, which can be used to determine where to limit the optimization depth.

  • Encoding/decoding schemes reduce the number of optimization layers but increase the circuit’s complexity by introducing more gates. Therefore, when selecting an encoding scheme, we must take into account the complexity of the generated circuit and the number of required gates, as well as the number of classical resources (memory, CPU) it will require. There must be a trade-off between circuit complexity and the desired problem accuracy.

  • Despite the fact that amplitude encoding provided the greatest accuracy, it could not be used to simulate a 12-qubit VRP scenario due to the large number of gates required. Angel encoding, on the other hand, was found to be much simpler due to a significantly smaller number of gates, as well as providing excellent accuracy (\(96\%\) for COBYLA, and \(92\%\) for SLSQP and L_BGFS_B in variable hamiltonian simulation) across all the available optimizers. This again demonstrates that the complexity of circuits and the number of gates used are the most important considerations when choosing an encoding/decoding scheme.

    Table 3 Table consisting of quantum depth and quantum cost of various encoding schemes compared with standard VQE implementation of VRP
    Full size table

  • It can be noticed that AgE performs the best in terms of circuit complexity and accuracy rates due to the formation of a single layer of superposition. In other encodings (HO, IqpE), we observe multi-layered complex superposition structures, which is the reason for fluctuations or error rates. Also in the fact that increasing layers also increases the superposition structures and therefore decreases the accuracy.

  • Using COBYLA as an optimizer, HO encoding yielded intriguing results with reduced accuracy in circuits with fewer qubits (6 qubits) and higher accuracy in circuits with more qubits (12 qubits) for both fixed and variable Hamiltonian simulations. The trend is disregarded by SLSQP and L_BGFS_B. This demonstrates that the algorithm’s performance is extremely dependent on the optimizer; therefore, when evaluating the algorithm’s performance, the most efficient optimizer should be selected by comparing the available optimizers.

  • The IQP encoding scheme performed the worst in this experiment, with the lowest accuracy and highest error rates among all other encodings used for one-layer, two-layer, fixed, and variable Hamiltonians simulations. Therefore, the IqpE method cannot be used to solve VRP using QSVM.

  • All of the optimizers used in the experiments performed well across AE, AgE, and HO encodings; however, COBYLA outperformed the other two due to its consistently high level of accuracy, but SLSQP is more resistant to accuracy fluctuations caused by an increase in optimization depth or in the presence of multi-layered circuits.

4.3 Complexity and cost considerations

Along with the inferences and observation, we want to touch upon the complexity of the quantum circuits particularly the number of qubits, quantum cost, and circuit depth in this section. Here, we compare the SVM VRP solution with the standard VQE solution that is described in our previous work (Mohanty et al. 2023).

In order to compare, we have listed the quantum cost and circuit depth of standard VQE implementation with that of SVM VRP implementation each consisting of 6 and 12 qubits with one and two layers (Table 3). We can broadly derive the following observations from Table 3.

  • Transitioning from the conventional VQE method for VRP to using a quantum support vector machine results in a rise in the number of gates, elevating quantum depth and cost.

  • The rise in quantum cost and depth is not linear and varies depending upon the encoding scheme employed when comparing different encoding schemes with the conventional VQE implementation of VRP.

  • Increasing the number of layers and qubits in various techniques consistently leads to a rise in quantum depth and quantum cost. The extent of this increase varies based on the encoding method used.

  • For the VQE standard, amplitude encoding, and angle encoding, there is a proportionate rise in depth and cost, roughly doubling as we go through layers and qubits.

  • Higher-order encoding (HO) and IQP encoding exhibit a much smaller rise in cost and depth over layers and qubits compared to amplitude and angle encoding.

  • As the number of qubits increases from 6 to 12, the depth and cost rise by almost four times compared to the 6-qubit scenario; however, the increase is somewhat less than four times for HO and IQP encodings.

  • Based on the depth and cost analysis, it is evident that HO and IQP provide lower costs and circuit depths compared to amplitude and angle encoding, but their findings are less precise. Amplitude and angle encoding provide the highest level of precision. This helps in determining the balance between cost and accuracy when choosing encoding techniques for QSVM.

4.4 Experimental setup, data gathering, and statistics

This experiment is conducted within the ambit of the QISKIT framework. While performing the experiment, we used a quantum instance object, and the ansatz runs inside the quantum instance object. A random seed is added to quantum instance to stabilize VQE results. All the experiments have been run \(50+50\) times, one with a fixed Hamiltonian matrix and the other by varying the Hamiltonian matrix. The objective of the experiments is to ensure that the results of experiments are just not dependent on a single Hamiltonian. This is also to ensure that the used circuits achieve classical minimum or near classical minimum regardless of the Hamiltonian used. Thus, apart from the plots, Tables 1 and 2 become the figure of merit. In addition to the many hours of testing and debugging, it is to be noted that the results reported here amounted to 150 h of CPU time on a 24-core AMD workstation using Qiskit’s built-in simulators (et al. 2024).

5 Conclusion

In this paper, we presented a novel technique for solving VRP through the use of a 6- and 12-qubit circuit-based quantum support vector machine (QSVM) with a variational quantum eigensolver for both fixed and variable Hamiltonians. In the experiment, multiple encoding strategies were used to convert the VRP formulation into a QSVM and solve it. In addition, we utilized multiple classical optimizers available within the QISKIT framework to measure the output variation and accuracy rates. Consequently, our machine learning-based approach to resolving VRP has proven fruitful thus far. Using a QSVM to implement a gate-based simulation of a 3-city or 4-city VRP on a 6-qubit or 12-qubit system accomplishes the goal. The method not only resolves VRP, but also outperforms the conventional method of resolving VRP via multiple optimization phases involving only VQE and QAOA. In addition, selecting appropriate encoding methods establishes the optimal balance between circuit complexity and optimization depth, thereby enabling multiple approaches to solve CO problems using machine learning techniques.