QUBO-Based Simulated Annealing Approach for the Shortest Path Problem with Applications to Urban Transportation Networks

1. Introduction

The shortest path problem has been extensively studied within graph theory. However, when we deal with dense graphs containing a large number of nodes and edges, the difficulty of the problem increases considerably. These types of structures appear frequently in large-scale networks, including applications such as route planning in transportation systems [1,2], routing in communication networks [3], robotic navigation [4,5], social network analysis [6,7], bioinformatics [8], among other areas [9,10,11,12]. In these practical applications, the complex network of connections between nodes generates additional computational challenges that intensify as the graph’s density increases.

This problem not only constitutes a fundamental component in various practical applications, but also offers a conceptual framework for addressing more complex problems in graph theory. For instance, the graph isomorphism problem [13] seeks to determine whether two graphs are structurally identical based on the paths between their nodes. It also proves essential in the minimum spanning tree (MST) problem [14], where the goal is to minimize the total weight of edges to connect all nodes, as well as in the maximum flow problem, which optimizes the flow of resources through a network from a source node to a sink node [15]. Additionally, in the traveling salesman problem (TSP), shortest path strategies help identify the most efficient route to visit all nodes [14]. Finally, path planning in robotics and autonomous systems, crucial for navigation in dynamic environments, relies heavily on these algorithms to calculate optimal routes.

The shortest path problem comes in four main variants, each tailored to specific use cases. The single-pair shortest path (SPSP) problem, which seeks the minimal route between two designated nodes, is the focus of this study. However, its related variants include the single-source shortest path (SSSP), single-destination shortest path (SDSP), and all-pairs shortest path (APSP), which together highlight the problem’s versatility. While APSP can be solved naively through iterative resolution of SPSP, specialized algorithms such as Dijkstra’s algorithm [16] or Floyd’s algorithm [17] frequently leverage structural symmetries to achieve polynomial-time efficiency for SSSP and APSP. Nevertheless, as graph density increases, these optimized algorithms face prohibitive computational costs, revealing fundamental limitations in classical approaches.

Traditional exact algorithms, such as Dijkstra, Bellman-Ford [18], and A* [18], guarantee optimality through greedy or dynamic programming strategies. While they perform efficiently on sparse graphs, their performance deteriorates in dense networks due to exponential growth in possible routes and high resource demands. This scalability issue, common to many NP-hard combinatorial optimization problems, motivates the adoption of heuristic and metaheuristic techniques that trade guaranteed optimality for tractability.

Metaheuristics such as genetic algorithms (GA) [19,20], ant colony optimization (ACO) [21,22], particle swarm optimization (PSO) [23,24], simulated annealing (SA) [25,26,27,28], and quantum annealing (QA) [29] have demonstrated strong potential by exploring large search spaces through nature-inspired mechanisms. Other approaches, including grey wolf optimization (GWO) [30,31], artificial bee colony (ABC) [32,33], and tabu search (TS) [34], further expand the toolkit for addressing NP-hard problems like the TSP [35]. However, in dense graphs, these methods risk premature convergence or high computational overhead, underscoring the need for frameworks that balance exploration with robust constraint handling.

Quadratic Unconstrained Binary Optimization (QUBO) provides such a framework by encoding nodes, edges, and constraints into a unified quadratic objective [36,37]. This formulation captures the entire solution space without iteratively enumerating paths, making it suitable for large-scale instances. When combined with metaheuristics like SA, which probabilistically explores the search space to avoid local minima, QUBO models can efficiently navigate high-dimensional landscapes, achieving a balance between exploration and exploitation. This integration, referred to as QUBO–SA, emerges as a promising methodology for solving the SPSP problem across diverse graph structures.

In this work, we evaluate the performance of the QUBO-SA framework for SPSP across synthetic graphs with controlled connectivity properties and a real-world urban transportation network. We provide both theoretical formulations and empirical evaluations to assess the solution quality, computational efficiency, and scalability characteristics of our approach through probabilistic performance metrics.

Summarizing, the main contributions of this article has been addressed as

• An evaluation of the performance of the QUBO-SA framework for SPSP across synthetic graphs with controlled connectivity properties and a real-world urban transportation network.
• The implementation of a case study in urban areas to make a simulation of the QUBO-SA approach using QUBO optimization.
• The theoretical formulations and empirical evaluations to assess the solution quality, computational efficiency, and scalability characteristics of our approach through probabilistic performance metrics.

The paper is organized as follows: Section 2 reviews existing literature on SPSP solutions. Section 3 introduces the QUBO model and SA algorithm. Section 4 presents the SPSP mathematical formulation and its QUBO-based reformulation. Section 5 reports the experimental results, comparative performance analysis, and potential extensions and future research directions. Section 6 summarizes the main findings and implications of this work.

2. Related Work

The Single-Pair Shortest Path (SPSP) problem, which seeks the shortest path between two specific nodes in a graph, is a special case of the more general Single-Source Shortest Path (SSSP) problem. Much of the research on shortest path algorithms focuses on SSSP due to its broader applicability, and many techniques for SSSP can be adapted to SPSP.

In unweighted graphs, Breadth-First Search (BFS) [38] provides an effective strategy, as it identifies the shortest route by minimizing the number of edges crossed. For weighted graphs, however, Dijkstra’s algorithm [16] remains a central technique, enabling the calculation of minimum path lengths from a chosen source to all other nodes. Originally $O\left(n^2\right)$, its complexity was improved to $O(m+nlogn)$ with advanced data structures such as Fibonacci heaps [39] and relaxed heaps [40], enabling faster insertions and key updates.

For directed graphs, specialized problems and optimizations have been explored. The bottleneck path problem, maximizing the minimum edge weight along a path, was addressed by Gabow and Tarjan [41] in $O(mlogn)$, later improved to $O\left(m\sqrt{logn}\right)$ by Duan et al. [42] via randomized methods. Williams [43] proposed an $O(mloglogn)$ algorithm for the non-decreasing path problem, and Duan et al. [44] recently broke the long-standing $O(m+nlogn)$ bound for SSSP on sparse directed graphs with non-negative weights, achieving $O\left(m{log}^{2/3}n\right)$ in the comparison–addition model.

Undirected graphs admit different optimizations. Pettie and Ramachandran [45] achieved $O\left(m\alpha \right(m,n)+min\{nlogn,nloglogr\left\}\right)$ in the comparison–addition model, where $\alpha$ is the inverse Ackermann function and r bounds edge weight ratios. For minimum spanning trees (MST), chih Yao [46] achieved $O(mloglogn)$, while Tarjan improved to $O\left(m\sqrt{logn}\right)$. These advances show how algorithmic strategies differ between directed and undirected graphs.

Quantum algorithms have also been applied to SPSP. Dürr et al. [47] proposed a quantized Dijkstra’s algorithm with complexity $O\left(\sqrt{n}m{log}^{2}n\right)$ using quantum minimum finding, and addressed the s–t connectivity problem with $\Theta \left(n\right)$ adjacency list queries. Building on this, Belovs and Reichardt [48] developed an $O\left(n\sqrt{L}\right)$ method for detecting an s–t path, where L is the longest path length, later improved by Belovs [49] to $O\left(\sqrt{RW}\right)$ queries, with R the effective resistance between s and t and W the total graph weight. Jeffery et al. [50] matched the Belovs–Reichardt detection time while reconstructing the full path in $\tilde{O}\left(L^{3/2}n\right)$ adjacency matrix queries. More recently, Wesołowski and Piddock [51] introduced bounded-error quantum algorithms in the adjacency list model for structured graphs with complexities $\tilde{O}\left(l^2\sqrt{m}\right)$ and $\tilde{O}\left(l\sqrt{m}\right)$, where l is the shortest-path length, and Ashvinkumar et al. [52] presented a reduction framework for handling negative edge weights in SSSP, achieving quantum edge-query complexity $m^{1/2}{n}^{1/2+o\left(1\right)}$ for directed graphs.

The approach proposed in this work addresses the Single-Pair Shortest Path (SPSP) problem in weighted undirected graphs through a Quadratic Unconstrained Binary Optimization (QUBO) formulation solved via Simulated Annealing. This method integrates the constraints of source, destination, and intermediate nodes into a unified quadratic cost function, ensuring that the only valid solutions minimize the cost while adhering to the shortest path constraints.

3. Preliminaries

3.1. Quadratic Unconstrained Binary Optimization (QUBO)

The Quadratic Unconstrained Binary Optimization (QUBO) problem is a fundamental framework in Combinatorial Optimization (CO) problems, involving the minimization of a quadratic objective function over binary variables. The problem is defined as follows:

$$f\left(x\right)=x^{T}Qx=\sum _{i,j=1}^n{Q}_{ij}{x}_i{x}_j,$$

(1)

where $x\in {\{0,1\}}^n$ represents a vector of n binary decision variables, and $Q\in {\mathbb{R}}^{n\times n}$ is a symmetric matrix. The elements $Q_{ij}$ capture the interaction between variables $x_i$ and $x_j$, with diagonal entries $Q_{ii}$ representing the individual weights of the variables, and off-diagonal entries $Q_{ij}$ indicating quadratic interactions between variable pairs. Due to the binary character of the variables, the product $x_i{x}_j$ is always either 0 or 1, and the symmetry of Q stems from the commutative property of multiplication ($Q_{ij}=Q_{ji}$). The objective function, being a quadratic polynomial, concentrates exclusively on second-degree terms, either involving products of distinct binary variables or squares of individual variables.

In its native form, QUBO is unconstrained, but constraints from real-world optimization problems can be embedded as penalty terms in the objective function. These penalties function to discourage invalid or infeasible solutions by raising the objective value when constraints are violated. Despite its flexibility in representing various optimization problems, the QUBO formulation is NP-hard, meaning that solving it optimally becomes computationally intractable for large problem instances [53,54].

3.2. Simulated Annealing for QUBO Models

Simulated Annealing (SA) represents a stochastic metaheuristic that was introduced by Kirkpatrick et al. [55], drawing inspiration from the physical annealing process in metallurgy. In this process, materials are heated and then gradually cooled to eliminate defects and achieve a stable crystalline structure with low energy. This physical concept translates naturally into the realm of combinatorial optimization, where ”energy” corresponds to the objective function value, and the search process seeks to identify a global optimum among numerous local minima. SA has proven successful across a wide range of combinatorial optimization (CO) problems, including Quadratic Unconstrained Binary Optimization (QUBO), where solution spaces expand exponentially with problem size and exhibit rugged, multi-modal landscapes.

The foundation of SA lies in the Metropolis acceptance criterion [56], which is a probabilistic rule derived from Monte Carlo methods in statistical mechanics. This criterion guides the iterative exploration of the solution space. Beginning with an initial solution $\mathbf{x}$, the algorithm generates a neighboring configuration ${\mathbf{x}}^{\prime }$ based on a predefined neighborhood structure. The objective function value difference (energy difference) between the candidate and current solutions is calculated as:

$$\Delta E=f\left({\mathbf{x}}^{\prime }\right)-f\left(\mathbf{x}\right),$$

(2)

where $f\left(\mathbf{x}\right)$ represents the QUBO objective function $f\left(\mathbf{x}\right)={\mathbf{x}}^{\top }Q\mathbf{x}$.

The decision to accept ${\mathbf{x}}^{\prime }$ is fundamentally based on the Boltzmann distribution, which describes the relative probability of a system occupying a state with energy E at temperature T:

$$P\left(E\right)\propto e^{-{\displaystyle \frac{E}{k_{B}T}}},$$

where $k_B$ denotes the Boltzmann constant. Within the optimization framework, the probability of accepting a move from $\mathbf{x}$ to ${\mathbf{x}}^{\prime }$ is given by:

$$P_{\mathrm{accept}}=e^{-{\displaystyle \frac{\Delta E}{T}}},$$

(3)

when $\Delta E>0$, whereas moves that improve the objective ($\Delta E<0$) are always accepted. The Boltzmann factor $exp(-\Delta E/T)$ ensures that when temperatures are high, the system regularly accepts inferior solutions, which promotes broad exploration of the search space. Conversely, as temperatures decrease, the system grows more selective and focuses on exploitation around promising regions.

Within the QUBO-SA framework, each candidate solution takes the form of a binary vector $\mathbf{x}\in {\{0,1\}}^n$ that encodes the problem’s decision variables. The quadratic coefficient matrix $Q\in {\mathbb{R}}^{n\times n}$ fully characterizes the energy landscape, and the objective function $f\left(\mathbf{x}\right)$ serves as the energy measure we aim to minimize. The complete procedure for applying SA to QUBO problems is outlined in Algorithm 1, which demonstrates how the main components work together systematically. The algorithm operates through the following steps:

1.

Initialization: The process begins by generating a random feasible solution $\mathbf{x}$, which acts as the starting point for exploration. This initial solution is also designated as the current best-known solution ${\mathbf{x}}^{*}$, providing a benchmark for tracking improvements. The initial temperature $T_0$ is established to regulate the probability of accepting uphill moves.

2.

Neighbor generation: A candidate solution ${\mathbf{x}}^{\prime }$ is produced by applying a neighborhood operator to $\mathbf{x}$. This typically involves flipping the state of one or more components in the binary vector, allowing the algorithm to examine nearby configurations.

3.

Acceptance criterion: The algorithm computes the energy difference $\Delta E=f\left({\mathbf{x}}^{\prime }\right)-f\left(\mathbf{x}\right)$. When $\Delta E<0$, the candidate solution ${\mathbf{x}}^{\prime }$ is immediately accepted. Otherwise, acceptance occurs with probability $exp(-\Delta E/T)$. If the solution is accepted and satisfies $f\left(\mathbf{x}\right)<f\left({\mathbf{x}}^{*}\right)$, then ${\mathbf{x}}^{*}$ is updated accordingly.

4.

Cooling schedule: Temperature T decreases progressively following a geometric schedule $T\leftarrow \alpha T$ where $\alpha \in (0,1)$, which maintains a balance between exploration and exploitation throughout the search. The procedure continues until temperature drops below a minimum threshold $T_{min}$ or reaches the maximum iteration count $k_{max}$, thereby ensuring the algorithm remains computationally feasible.

The QUBO-SA algorithm integrates these components into a unified iterative process. At each step, a candidate solution ${\mathbf{x}}^{\prime }$ is generated from the current state $\mathbf{x}$ and evaluated with respect to the QUBO objective. The acceptance rule balances diversification and intensification by allowing occasional uphill transitions, while improvements are tracked through the best-known solution ${\mathbf{x}}^{*}$. As the temperature parameter decreases, the search progressively shifts from broad exploration of the energy landscape toward focused refinement around promising regions. This combination of probabilistic exploration and controlled convergence defines the solver developed in this work.

4. Methods

This section describes the mathematical framework that underlies the proposed approach. We start by defining the Single-Pair Shortest Path problem in undirected weighted graphs, establishing the notation and structural constraints that define valid solutions. Following this, we adopt the QUBO formulation for the shortest path problem proposed by Krauss and McCollum [57], which consolidates both the path cost minimization objective and the structural validity constraints into a unified quadratic function over binary variables. Finally, we detail the algorithmic construction of the quadratic coefficient matrix Q, which enables the application of Simulated Annealing for efficient exploration of the solution space across diverse graph topologies.

4.1. Single-Pair Shortest Path Problem (SPSP)

Let $G=(V,E,C)$ represent an undirected graph, where $V=\{1,2,\dots ,n\}$ is a finite set of vertices, $E\subseteq \left\{\right(i,j)\mid i,j\in V,i\ne j\}$ is a set of edges, and $C=\{c_{i,j}\mid (i,j)\in E\}\subseteq \mathbb{R}$ is a set of real-valued edge costs associated with each edge $(i,j)$. In this article, we assume that all edge costs are non-negative, i.e., $c_{i,j}\ge 0$ for all $(i,j)\in E$, and therefore, negative edge costs are not considered.

In an undirected graph, the edges do not have a specific direction, meaning that if $(i,j)$ is an edge in the graph, then $(j,i)$ is also an edge. The connection between the vertices i and j is bidirectional.

A path in the graph from a node s to a node t is defined as a finite sequence of edges $(e_1,e_2,\dots ,e_k)$ where each edge $e_i=(v_{i-1},v_i)$ are adjacent. The shortest path problem involves finding such a sequence of edges starting at a source node s and ending at a target node t, where the starting node of the first edge in the sequence is the source point s, the arrival node of the last edge in the sequence is the target point t, and for each edge in the sequence, the arrival node of the previous edge is the same as the starting node of the next edge.

A path (P) from the source node s to the target node t can be expressed as:

$$\mathrm{P}=\{(v_0,v_1),(v_1,v_2),\dots ,(v_{k-1},v_k)\}$$

(4)

where $v_0=s$ and $v_k=t$. The total cost $C\left(P\right)$ of P is determined by the sum of the weights of its constituent edges:

$$C\left(\mathrm{P}\right)=\sum _{i=1}^k{c}_{v_{i-1},v_i}$$

(5)

The objective of the shortest path problem is to find the path that minimizes the total cost $C\left(\mathrm{P}\right)$.

It is important to note that the shortest path between two nodes is not necessarily unique. Multiple distinct paths ${\mathrm{P}}_1,{\mathrm{P}}_2,\dots ,{\mathrm{P}}_m$ may exist between the same pair of nodes such that $C\left({\mathrm{P}}_1\right)=C\left({\mathrm{P}}_2\right)=\dots =C\left({\mathrm{P}}_m\right)$, all of which would be considered optimal solutions satisfying $C\left({\mathrm{P}}_i\right)=C\left({\mathrm{P}}^{*}\right)$. For example, as illustrated in Figure 1, consider a scenario where the graph contains two alternative paths: ${\mathrm{P}}_1=\{(x_s,x_3),(x_3,x_t)\}$ and ${\mathrm{P}}_2=\{(x_s,x_4),(x_4,x_t)\}$. If both paths yield $C\left({\mathrm{P}}_1\right)=C\left({\mathrm{P}}_2\right)=11$, then either constitutes a valid optimal solution.

The complexity of solving the shortest path problem can vary significantly depending on the structure of the graph and the nature of the edge costs. Classical algorithms deliver efficient polynomial-time solutions for numerous cases; however, problem complexity can escalate dramatically under more demanding scenarios, especially when dealing with additional constraints or large-scale graphs. Under such circumstances, advanced optimization techniques become necessary to overcome the computational challenges that arise.

4.2. QUBO Model Formulation

The SPSP problem can be expressed as a Quadratic Unconstrained Binary Optimization (QUBO) model, where both the path cost objective and structural feasibility requirements are jointly represented through a quadratic energy function defined over binary decision variables. This encoding approach utilizes linear coefficients for representing edge weights and quadratic penalty coefficients for enforcing topological constraints on the solution structure.

A valid path configuration needs to satisfy three essential requirements: (i) it must originate at the designated source node s, (ii) it must terminate at the specified target node t, and (iii) it must maintain edge connectivity throughout the entire path, thereby preventing fragmentation or branching at any vertex. To encode these conditions, the formulation uses $\left|V\right|+\left|E\right|-2$ binary variables, which are divided into $\left|V\right|-2$ node indicators representing intermediate vertices and $\left|E\right|$ edge indicators covering all graph connections.

The variable assignment follows a structured convention: each intermediate node $i\in V\setminus \{s,t\}$ corresponds to a binary indicator $x_i$, while each edge $(i,j)\in E$ is assigned a binary selector $x_{ij}$ denoting its inclusion in the path. The vector therefore represents a complete system state:

$$\mathbf{x}=(x_1,\dots ,x_{\left|V\right|-2},\dots ,x_{i,j},\dots ),$$

where the source and target nodes are implicitly activated in all configurations. The objective function assumes the canonical QUBO form:

$$Q\left(\mathbf{x}\right)=Q_c\left(\mathbf{x}\right)+\lambda Q_p\left(\mathbf{x}\right),$$

(6)

The cost component $Q_c\left(\mathbf{x}\right)$ aggregates edge weights across the selected path:

$$Q_c\left(\mathbf{x}\right)=\sum _{i\in V\setminus \{s,t\}}{c}_i{x}_i+\sum _{(i,j)\in E}{c}_{i,j}{x}_{i,j},$$

(7)

while the constraint component $Q_p\left(\mathbf{x}\right)$ enforces structural validity through three penalty mechanisms:

$$Q_p\left(\mathbf{x}\right)=Q_s\left(\mathbf{x}\right)+Q_t\left(\mathbf{x}\right)+Q_{\mathrm{cont}}\left(\mathbf{x}\right).$$

(8)

The penalty factor $\lambda$ serves a critical function in ensuring that constraint violations dominate the objective value, effectively preventing infeasible solutions from being considered optimal. A properly calibrated penalty factor needs to satisfy $\lambda >c_P$, where $c_P$ denotes the cost of the shortest path. Given that $c_P$ remains unknown beforehand, we employ a conservative upper bound: $\lambda ={\sum }_{(i,j)\in E}\left|c_{i,j}\right|$, representing the sum of all edge costs in absolute value. This selection ensures that any configuration violating structural constraints receives a penalty that exceeds the maximum possible path cost, thereby effectively eliminating such configurations from consideration as optimal solutions. The scaling mechanism guarantees that the global minimum of $Q\left(\mathbf{x}\right)$ corresponds exclusively to valid path structures exhibiting minimal cost.

The source constraint $Q_s\left(\mathbf{x}\right)$ guarantees that precisely one edge is incident to the source node s:

$$Q_s\left(\mathbf{x}\right)=-x_s^2+{\left(x_s-\sum _{j\in \delta \left(s\right)}{x}_{s,j}\right)}^2,$$

(9)

where $\delta \left(s\right)$ represents the neighborhood of s and $x_s=1$ by construction. This penalty function imposes the requirement that a valid path must begin from exactly one outgoing edge at the source vertex. Upon expansion, the quadratic term produces:

$${\left(x_s-\sum _{j\in \delta \left(s\right)}{x}_{s,j}\right)}^2=x_s^2-2x_s\sum _{j\in \delta \left(s\right)}{x}_{s,j}+{\left(\sum _{j\in \delta \left(s\right)}{x}_{s,j}\right)}^2.$$

(10)

The leading negative term $-x_s^2$ functions to compensate for the constant contribution when the constraint is met, yielding a global minimum of $-1$ achieved precisely when ${\sum }_{j\in \delta \left(s\right)}{x}_{s,j}=1$. Any departure from this condition, whether through zero edges (disconnected path) or multiple edges (branching), raises the penalty value, thus discouraging such configurations during optimization. The quadratic structure guarantees that violations receive penalties proportional to their magnitude, with the squared difference amplifying deviations from the required degree.

The target constraint $Q_t\left(\mathbf{x}\right)$ enforces a symmetric condition at the destination node t:

$$Q_t\left(\mathbf{x}\right)=-x_t^2+{\left(x_t-\sum _{i\in \delta \left(t\right)}{x}_{i,t}\right)}^2,$$

(11)

with expansion

$${\left(x_t-\sum _{i\in \delta \left(t\right)}{x}_{i,t}\right)}^2=x_t^2-2x_t\sum _{i\in \delta \left(t\right)}{x}_{i,t}+{\left(\sum _{i\in \delta \left(t\right)}{x}_{i,t}\right)}^2,$$

(12)

where $x_t=1$ and $\delta \left(t\right)$ denotes the neighborhood of t. This constraint ensures that the path concludes through exactly one incoming edge at the target vertex, mirroring the functionality of the source constraint. The minimum value of $-1$ is reached when ${\sum }_{i\in \delta \left(t\right)}{x}_{i,t}=1$, thereby ensuring proper path termination. Similar to the source constraint, the negative quadratic offset modifies the baseline energy so that valid configurations produce negative contributions, whereas infeasible terminations, whether through no incoming edges or multiple simultaneous entries, incur positive penalties that render them ineligible for optimality.

The path continuity constraint $Q_{\mathrm{cont}}\left(\mathbf{x}\right)$ guarantees that each intermediate vertex maintains connectivity through exactly two incident edges:

$$Q_{\mathrm{cont}}\left(\mathbf{x}\right)=\sum _{i\in V\setminus \{s,t\}}{Q}_i\left(\mathbf{x}\right),$$

(13)

where each nodal contribution is expressed as

$$Q_i\left(\mathbf{x}\right)={\left(2x_i-\sum _{j\in \delta \left(i\right)}{x}_{i,j}\right)}^2,$$

(14)

with explicit expansion

$${\left(2x_i-\sum _{j\in \delta \left(i\right)}{x}_{i,j}\right)}^2=4x_i^2-4x_i\sum _{j\in \delta \left(i\right)}{x}_{i,j}+{\left(\sum _{j\in \delta \left(i\right)}{x}_{i,j}\right)}^2.$$

(15)

Here, $\delta \left(i\right)$ encompasses all edges incident to node i. This constraint implements the essential principle that intermediate nodes on a simple path must display degree two: one edge for entrance and one for exit. When $x_i=0$, the node is excluded from the path, and the penalty drops to zero provided no incident edges are selected (${\sum }_{j\in \delta \left(i\right)}{x}_{i,j}=0$). In contrast, when $x_i=1$, signaling inclusion in the path, the penalty is minimized at zero only when ${\sum }_{j\in \delta \left(i\right)}{x}_{i,j}=2$. Any violation, such as a degree-one node creating a dead end, or a degree-three (or higher) node causing branching, produces a positive penalty proportional to the squared deviation from the required degree. The summation across all intermediate vertices ensures that path continuity is enforced globally throughout the entire configuration. Taken together, the minimization of $Q_p\left(\mathbf{x}\right)$ ensures that optimal configurations correspond to structurally valid paths connecting s to t without internal fragmentation or extraneous branches.

4.3. Algorithm Description

The construction of the quadratic coefficient matrix Q for the SPSP problem is achieved through a systematic process distributed across two algorithms. Algorithm 2 establishes the variable encoding and edge cost representation ($Q_c$), whereas Algorithm 3 enforces path validity through quadratic penalty terms corresponding to source and target nodes ($Q_s$ and $Q_t$), intermediate nodes ($Q_i$), and general path continuity ($Q_{\mathrm{cont}}$).

The procedure in Algorithm 2 begins by creating a bijective mapping ${\varphi }_E$ that associates each edge $(i,j)\in E$ with a unique binary variable index. Given that the graph is undirected, both orientations $(i,j)$ and $(j,i)$ map to the same index. Simultaneously, for each node $v\in V$, the set $\delta \left(v\right)$ of incident edge indices is compiled, enabling direct access to all edges connected to that vertex. This establishes the variables corresponding to nodes and edges, which will subsequently interact through the constraints. The variable vector incorporates binary variables for all vertices, including the source and destination nodes; later, $x_s=x_t=1$ are fixed through the penalty terms that enforce their activation within the feasible path.

Following index assignment, the edge costs $c_{i,j}$ are positioned on the diagonal entries $Q[k+n,k+n]$, where $k={\varphi }_E\left[(i,j)\right]$. This encodes the path cost directly into the matrix as $Q_c$, ensuring that the total path length appears as the sum of selected edge variables.

To guarantee structural validity of the path, constraints are implemented through Algorithm 3. Each node is handled according to its target degree: $d_v=1$ for the source node s and the target node t, and $d_v=2$ for intermediate nodes. These operations modify Q in-place, introducing quadratic and pairwise terms that realize the structural penalties included in $Q_p$. The resulting matrix Q now incorporates both the edge costs ($Q_c$) and the node and edge penalty terms ($Q_p$).

Algorithm 3 implements the constraint enforcement mechanism. For a node v with target degree $d_v$, three types of penalties are applied: (i) a quadratic penalty $d_v^2\lambda$ is added to $Q[v,v]$ (line 2), representing the cost of selecting the node itself; (ii) for each edge $e\in \delta \left(v\right)$, the diagonal entry $Q[e+n,e+n]$ is incremented by $\lambda$, and the coupling $Q[v,e+n]$ is decremented by $2d_v\lambda$ (lines 4–5), encouraging configurations that satisfy the node degree; (iii) for each pair of edges $(e,f)\in \delta \left(v\right)$ with $e<f$, a pairwise term $2\lambda$ is added to $Q[e+n,f+n]$ (lines 7–9), penalizing configurations where multiple edges would violate the target degree. These penalties correspond to the components $Q_{\mathrm{cont}}$ and $Q_i$, ensuring continuity and correct degree enforcement for each node.

5. Experimental Process and Results

This section describes the experimental framework and results obtained from evaluating the proposed QUBO-SA algorithm in comparison with the classical Dijkstra solver [16].

The comparison focuses on two main aspects: (i) solution reliability, captured through the empirical success probability $p_{\mathrm{success}}$, and (ii) computational efficiency, assessed using the Time-to-Solution ($\mathrm{TTS}$) metric and the relative runtime ratio $R\left(p_{\mathrm{target}}\right)$.

Computational Resources. All experiments were performed using Julia 1.11.5 on hardware equipped with an AMD Ryzen AI 9 HX 370 processor (12 cores / 24 threads, base clock 2.0 GHz) and 32 GB of RAM. No GPU acceleration was utilized.

5.1. Experimental Design and Evaluation Metrics

The evaluation procedure consists of three primary stages: graph preparation, QUBO encoding and optimization via simulated annealing (SA), and performance assessment based on probabilistic reliability. Both synthetic and real-world graph instances were employed to investigate how graph size and density affect solver behavior.

Graph preparation: Two categories of graph instances were considered:

• Synthetic graphs: Undirected weighted graphs were constructed with a target average degree k, from which the corresponding density d is calculated. This parametrization facilitates systematic exploration of sparse ($k=1$), moderately dense ($k=2$), and dense ($k=3$) connectivity regimes across multiple node counts ($n=10,20,30,40$). The complete generation procedure and pseudocode are detailed in Appendix A.
• Urban Transportation Network. 
  This real-world urban transportation network was derived from the downtown area of Querétaro City, México, a densely populated district with intricate routing demands arising from the coexistence of residential, commercial, and administrative zones. In this graph, nodes represent intersections or decision points, whereas edges denote bidirectional road segments weighted by geodesic distances calculated via the Haversine formula. The resulting network encapsulates realistic spatial constraints, including irregular connectivity patterns and hierarchical road organization characteristic of urban infrastructures.
  Figure 2 offers a representation of the study area: subfigure 2a positions Querétaro within its broader geographic context through a progressive zoom from national to municipal scale, whereas subfigure 2b illustrates the extracted network topology of the downtown area. In the latter, node coordinates match actual intersection positions obtained from OpenStreetMap data, and the designated source (green) and destination (blue) nodes were chosen to define a routing scenario that traverses multiple urban blocks, thereby encouraging path diversity across the sparse infrastructure.
  The extracted network contains $n=443$ nodes and $m=245$ edges, yielding an average degree of $k\approx 1.11$ and a density of $d\approx 0.25\%$. This structural sparsity is typical of urban road systems, where intersections generally connect to only one or two neighboring streets, forming quasi-tree topologies. Such configurations impose a significant topological constraint on path diversity, thereby testing the capacity of the QUBO formulation to preserve connectivity continuity under limited redundancy.
  The routing scenario established for this network involves navigating 37 edges with an optimal path cost of 1,649 meters, representing the shortest valid route between the designated source and destination nodes. Table 1 summarizes the main topological and spatial characteristics of the Querétaro network that served as the real-world benchmark throughout the experimental evaluation.

QUBO-SA encoding and solver configuration 

Each graph instance was transformed into a QUBO matrix following the formulation described in Section 4.2, where both the objective and constraint terms are incorporated within the quadratic cost coefficients. The optimization process adheres to the simulated annealing (SA) approach, as outlined in Algorithm 1, which utilizes a probabilistic acceptance of new states and a geometric cooling schedule. The solver parameters were calibrated as indicated in Table 2, guaranteeing consistent performance across different graph sizes.

Evaluation metrics 

The performance assessment of the QUBO-SA solver depends on probabilistic indicators that capture both solution reliability and computational efficiency, recognizing the inherently stochastic character of simulated annealing. Instead of depending on single-run runtimes or best-case outcomes, we embrace a statistical framework where repeated executions yield an empirical success probability $p_{\mathrm{success}}$, defined as the proportion of runs that successfully recover the optimal path within the permitted iteration budget. This probability measures solver robustness and consistency across multiple trials, providing a more informative metric than isolated runtime values.

To relate this probabilistic measure to expected computational effort, we utilize the Time-to-Solution ($\mathrm{TTS}$) metric [58,59,60,61], defined as:

$$\mathrm{TTS}\left(p_{\mathrm{target}}\right)=t_{\mathrm{run}}·{\displaystyle \frac{ln(1-p_{\mathrm{target}})}{ln(1-p_{\mathrm{success}})}},$$

where $t_{\mathrm{run}}$ represents the average runtime per single QUBO-SA execution, and $p_{\mathrm{target}}$ indicates the desired overall confidence of achieving an optimal solution at least once. Intuitively, $\mathrm{TTS}$ estimates the expected total time necessary to reach a prescribed reliability level by considering the number of repetitions required to attain that success probability. A solver with high $p_{\mathrm{success}}$ naturally produces smaller $\mathrm{TTS}$, whereas one with erratic convergence (low $p_{\mathrm{success}}$) will display exponentially increasing $\mathrm{TTS}$ values even when individual runs are fast.

In this study, we chose $p_{\mathrm{target}}=0.99$ to represent high-confidence convergence aligned with standard benchmarks in stochastic optimization literature. Furthermore, we report values for $p_{\mathrm{target}}=0.90$ to reveal the trade-off between reliability and computational efficiency: lowering the target confidence typically reduces $\mathrm{TTS}$, indicating that slightly less stringent reliability requirements can result in substantial runtime savings.

For comparative evaluation, the deterministic Dijkstra algorithm functions as the baseline reference, supplying a mean runtime $t_{\mathrm{D}}$ calculated over 100 independent executions. This enables the definition of the relative runtime ratio:

$$R\left(p_{\mathrm{target}}\right)={\displaystyle \frac{\mathrm{TTS}\left(p_{\mathrm{target}}\right)}{t_{\mathrm{D}}}},$$

which expresses the time required by the probabilistic solver to reach a given confidence level relative to the runtime of the deterministic benchmark. By normalizing both runtimes, R provides a scale-invariant comparison across heterogeneous graph instances and hardware environments: $R<1$ indicates that QUBO-SA attains the desired reliability faster than Dijkstra, whereas $R>1$ implies additional computational effort to achieve equivalent confidence.

5.2. Comparative Performance Across Graph Regimes

Table 3 summarizes the relative runtime ratios $R_{0.99}$ and $R_{0.90}$ across synthetic and urban networks. Confidence intervals ($R_{0.99}^{\mathrm{CI}}$) were obtained by propagating the binomial uncertainty of $p_{\mathrm{success}}$ through the $\mathrm{TTS}$ definition.

The synthetic benchmarks reveal a clear progression from near-trivial reliability at small scales to increasingly costly convergence as n grows. For $n=10$, the solver attains $p_{\mathrm{success}}\approx 1$ in two of three configurations, collapsing $\mathrm{TTS}$ to zero and yielding $R=0$; even in the less favorable case ($k=3$, $d=0.333$), $\mathrm{TTS}$ remains at 0.5 $\mu s$ with modest overhead relative to Dijkstra. By $n=20$, success probabilities drop into the 0.82–0.90 range and $\mathrm{TTS}$ increases to values between 65 and 185 $\mu s$, pushing $R_{0.99}$ into 2.09–2.82. This trend intensifies at $n=30$ and $n=40$, where $p_{\mathrm{success}}$ drops to 0.53–0.74 and $\mathrm{TTS}$ increases to values between 1,310 and 4,950 $\mu s$, producing $R_{0.99}=3.55$–$5.86$ with expanding confidence intervals, demonstrating that large instances amplify stochastic variability and the number of repetitions required for high-confidence convergence. Across all sizes, lowering the target from $p_{\mathrm{target}}=0.99$ to 0.90 consistently decreases runtime by approximately a factor of two ($R_{0.90}\approx \frac{1}{2}{R}_{0.99}$), suggesting that the primary cost stems from enhancing reliability rather than from single-run execution time. Taken together, the scaling with n and the associated decrease in d for fixed k constrain feasible path redundancy and hinder annealing transitions, shifting synthetic graphs from a regime of effectively deterministic recovery at small n to one where robustness requires exponentially more repetitions as n increases.

For the Querétaro urban transportation network ($n=443$, $d=0.0025$), the solver reached a success probability of $p_{\mathrm{success}}=0.57$, which translates into a $\mathrm{TTS}$ of 9,250 $\mu s$ compared to 6,850 $\mu s$ for Dijkstra. The resulting runtime ratio of $R_{0.99}=1.35$ shows that QUBO-SA requires about one third more time than the deterministic baseline to guarantee 99% confidence. When the confidence level is reduced to $p_{\mathrm{target}}=0.90$, the ratio decreases to $R_{0.90}=1.08$, revealing that most of the additional cost is focused on achieving very high reliability rather than on generating an optimal solution per se. This behavior aligns with the synthetic benchmarks, where the gap between $R_{0.99}$ and $R_{0.90}$ also captures the trade-off between robustness and efficiency.

Figure 3 displays the optimal 37-edge path of total cost 1,649 m, linking the source and target nodes across the urban district.

The structural characteristics of the urban graph help clarify these results. Its low density ($d=0.0025$) diminishes the number of feasible paths, which simplifies constraint satisfaction but also restricts redundancy, making it more challenging for the annealing process to escape local minima. Nevertheless, the observed runtime ratios remain competitive. Despite the substantially larger scale ($n=443$), the solver’s performance is comparable to that of synthetic graphs with $n=20$, where $R_{0.99}$ ranges from $2.09$ to $2.82$. This suggests that the spatial organization of real road networks, characterized by clustering around intersections and hierarchical arterial structures, introduces regularities that partially compensate for the computational burden of size, enabling QUBO-SA to navigate the solution space more efficiently than in purely random topologies.

5.3. Further Works

The experimental results described in this study reveal several promising avenues for future research. The neighbor generation mechanism in the current QUBO-SA implementation employs a simple bit-flip strategy, which may limit the efficiency of exploration in large solution spaces. An adaptive neighbor generation scheme that dynamically modifies the perturbation intensity based on search progress could substantially improve convergence rates. This mechanism could incorporate multiple neighborhood structures, choosing among them based on the current temperature regime or recent acceptance history, thereby strengthening the balance between diversification and intensification throughout the annealing schedule.

The penalty factor $\lambda$ is essential for guaranteeing constraint satisfaction. However, its current definition as the sum of all edge costs functions as a conservative upper bound that may introduce unnecessary numerical scaling in heterogeneous graphs. Performing a systematic sensitivity analysis to examine the relationship between $\lambda$, graph structure, and solution quality would yield valuable insights into optimal penalty calibration strategies. This analysis should investigate how different penalty scaling approaches influence both feasibility enforcement and convergence speed across graphs with varying cost distributions.

The QUBO formulation naturally extends to quantum annealing platforms, where the quadratic structure can be directly mapped onto quantum processing units, such as those manufactured by D-Wave Systems. Empirical evaluation on quantum hardware would permit a direct comparison of thermal and quantum annealing dynamics for the Shortest Path Selection Problem (SPSP), potentially uncovering quantum advantages in specific graph regimes or problem scales. Such an investigation would need to address hardware connectivity constraints, embedding overhead carefully, and the noise characteristics inherent to current quantum annealers.

6. Conclusions

This work has assessed the performance of a QUBO-based Simulated Annealing approach for solving the Single-Pair Shortest Path problem across diverse graph topologies. The experimental evaluation included synthetic graphs with controlled connectivity regimes spanning from sparse to dense configurations, alongside a real-world urban transportation network. Results indicate that the QUBO-SA framework achieves competitive performance relative to classical deterministic algorithms when reliability requirements are adjusted to practical confidence levels. For small synthetic instances, the solver reached near-perfect success rates with minimal overhead, whereas larger graphs displayed reduced success probabilities requiring additional computational effort to ensure high-confidence convergence. The urban network evaluation showed that structural characteristics of real-world graphs, including spatial clustering and hierarchical organization, partially offset the computational burden associated with problem scale, producing performance comparable to moderately sized synthetic instances despite significantly larger node counts.

The probabilistic nature of the solver introduces a fundamental trade-off between solution reliability and computational efficiency, as measured through the Time-to-Solution metric. Reducing the target confidence from 99% to 90% consistently lowered runtime overhead by approximately half across all tested configurations, demonstrating that most computational cost originates from achieving very high reliability rather than from the intrinsic difficulty of finding optimal solutions. This characteristic renders the approach particularly suitable for applications where occasional suboptimal solutions are acceptable or where multiple independent runs can be executed in parallel. Future work exploring adaptive neighborhood strategies, penalty calibration techniques, quantum hardware implementations, and alternative metaheuristic solvers will further clarify the practical boundaries of this approach and identify the specific problem regimes where QUBO-based optimization offers genuine advantages over traditional methods.