COMBINATORIAL OPTIMIZATION WITH D-WAVE OCEAN SDK
Solving Traveling Salesman Problems with D-Wave Quantum Annealing
Visit Every City For the Most Efficient Sales Trip!
Backstory
Imagine you’re a salesperson tasked with visiting several cities and returning to your starting point. You want to take the shortest route possible to save time and fuel. Here’s the catch: you must visit each city exactly once.
This is the Traveling Salesman Problem (TSP): finding the most efficient route that visits all destinations and returns home. The complexity scales factorially (n!) as the number of nodes increases, making brute-force methods infeasible for large instances and challenging even the fastest computers to solve efficiently. TSP is classified as an NP-hard problem because no polynomial-time algorithm is known for finding the optimal solution.
Example
We have cities (vertices, nodes): V = {A, B, C, D} and routes (edges) with weights:
- (A, B) = 10
- (A, C) = 15
- (A, D) = 20
- (B, C) = 35
- (B, D) = 25
- (C, D) = 30
We want to visit all cities exactly once and return to our starting point.
Suppose we start from city A: the shortest route is A → B → D → C and back to A. This is an example of a closed-loop TSP with a specified starting city. We could start with other cities and still have the same total cost, with different orders of the cities in the output.
For open-path TSP, the starting and potentially ending cities are given as part of the problem, which would affect the total cost because there is no return leg to the starting city.
In either case, we should ensure that no smaller loops exist that do not visit all cities. The below example shows a loop that only visits some cities.
This kind of loop is called a “subtour”.
Mathematical Definition
The mathematical formulation of an open-path TSP without a starting city condition can be written as an optimization problem:
Configuration:
- A set of cities: V = {1, 2, …, n}
- A distance matrix d[i, j], where d[i, j] is the distance (or cost) between city i and city j.
Variables:
x[i, j] is a binary variable, where:
Objective Function:
Minimize the total distance (or cost):
The total travel distance is the sum of the distances d[i, j] for all pairs of cities i and j, weighted by x[i, j], which equals 1 only if the route travels from city i to city j. By enforcing j≠i, we exclude self-loops from the optimization.
Since minimizing the sum has a trivial solution where all x[i, j] = 0, resulting in a zero-cost solution, we must add constraints to enforce a valid solution.
Constraints:
We ensure each city is visited exactly once (entered once and exited once):
Constraint 1: Each City is Entered Exactly Once
This ensures every city j has exactly one incoming route (i into j).
Constraint 2: Each City is Exited Exactly Once
This ensures every city i has exactly one outgoing route (j out of i).
Constraint 3: Subtour Elimination
The solution must form a single loop visiting all cities. In other words, we disallow smaller loops that do not visit all cities. Therefore, additional constraints are required to prevent disjoint loops (subtours).
There are two common approaches:
Approach 1: Set-Based Subtour Elimination (Dantzig–Fulkerson–Johnson formulation)
This constraint ensures that no small group (subset S) of cities forms its own loop by limiting the number of edges within that group to at most (|S| — 1). In other words, for any subset of cities S (with at least 2 cities but smaller than the total number of cities), there cannot be a way to connect all the cities in S without leaving the group.
For example, if 3 cities have 3 edges connecting them, they form a loop, which this constraint prevents.
In short, this constraint stops smaller groups of cities from forming disconnected loops, ensuring the entire set of cities stays connected in one continuous tour.
However, this form of subtour elimination is impractical for large-scale problems because the number of subsets grows exponentially with the number of cities:
For example, for n = 10, there are 1,012 subsets to consider.
import math
# Number of cities
n = 10
# Calculate the number of subsets of size 2 to n-1
num_subsets = 0
for i in range(2, n):
num_subsets += math.comb(n, i)
print(f"Number of subsets of size {size}: {num_subsets}") # prints 1012For n = 20, it’s 1,048,554. You get the idea.
Approach 2: Flow-Based Subtour Elimination (Miller-Tucker-Zemlin formulation)
where u[i] for i ∈ {2, 3, …, n} is an auxiliary variable representing the “position” of city i in the tour. These variables are introduced to track the order in which cities are visited:
- u[i] = 0 for the starting city.
- u[i] = 1 for the first city visited after the starting city.
- u[i] = 2 for the second city visited, and so on…
- u[i] = n − 1 for the last city visited before returning to the starting city.
When x[i, j] = 1 (meaning there is a direct route from city i to city j), then the constraint becomes:
This ensures that city j appears later in the tour than city i.
When x[i, j] = 0, the terms n⋅x[i, j] vanishes, and the constraint is trivially satisfied:
The flow constraints ensure that u[i] values are consistent with the travel order, eliminating disconnected subtours.
Note: This constraint won’t work for closed-loop TSPs because returning to the starting city violates the ordering condition. To resolve this, we fix a starting city and exclude the return to the start city from the constraint.
Environment Setup
We’ll use D-Wave’s Ocean SDK to solve the combinatorial optimization.
Install Required Libraries
I’d recommend setting up a local environment for this project using venv or conda first. The below is for venv:
# create a project folder
mkdir ocean_test
cd ocean_test# setup and activate a virtual environment
python3 -m venv venv
source venv/bin/activate# upgrade to obtain the latest libraries
pip install --upgrade pipThen, install the following Python libraries:
pip install dimod dwave-ocean-sdk matplotlib pandasI’ll use Jupyter notebook (.ipynb) to write the following program.
Defining the Problem
We define the problem with edge weights:
# Example TSP setup
# cities = [city1, city2, ...]
cities = ['A', 'B', 'C', 'D']
num_cities = len(cities)
# distances = {(city1, city2): distance}
distances = {
(0, 1): 10,
(0, 2): 15,
(0, 3): 20,
(1, 2): 35,
(1, 3): 25,
(2, 3): 30
}Solution 1: D-Wave NetworkX
We can solve TSP problems using D-Wave NetworkX, which is an extension of NetworkX.
First, we create a graph using NetworkX:
import networkx as nx
# Create a graph
G = nx.Graph()
for (i, j), d in distances.items():
G.add_edge(i, j, weight=d)Note: nx.Graph is undirected. So, (i, j) and (j, i) are the same edge.
Then, we can visualize it:
import networkx as nx
# Create a graph
G = nx.Graph()
for (i, j), d in distances.items():
G.add_edge(i, j, weight=d)
# Set labels for cities and edges
city_labels = {i: city for i, city in enumerate(cities)}
edge_labels = nx.get_edge_attributes(G, 'weight')
# Draw the graph with spring layout
spring_positions = nx.spring_layout(G, weight=None, seed=42)
import matplotlib.pyplot as plt
# Draw the graph with spring layout adjusted by edge weights
plt.figure(figsize=(8, 8))
nx.draw(G, spring_positions, with_labels=True, labels=city_labels,
node_size=700, node_color="lightblue", font_size=14)
nx.draw_networkx_edge_labels(G, spring_positions, edge_labels,
font_color='red', font_size=12)
plt.title("Traveling Salesman Problem Graph (Edge Length Adjusted by Weight)",
fontsize=16)
plt.show()
Let’s solve this TSP problem using ExactSolver, which will try all permutations to find the best solution (suitable for small problems like this one).
import dimod
import dwave_networkx as dnx
# Use ExactSolver to check all permutations
sampler = dimod.ExactSolver()
# Solve it!
answer = dnx.traveling_salesman(G, sampler)
print(answer)This returns [2, 0, 1, 3] (C → A → B → D) as the best solution, which is correct.
We can also specify the start node as A (=0):
answer = dnx.traveling_salesman(G, sampler, start=0)
print(answer)This returns [0, 1, 3, 2] (A → B → D → C), which is again correct, and the route has a different sequence of cities starting with A (=0).
Solution 2: Constrained Quadratic Model (CQM)
Here, we’ll use CQM to solve the same TSP.
Setting Up the Symmetric Distances
The original graph is meant for symmetric in terms of distances but it is set up like a DAG (directed acyclic graph).
# Example TSP setup
# cities = [city1, city2, ...]
cities = ['A', 'B', 'C', 'D']
num_cities = len(cities)
# distances = {(city1, city2): distance}
distances = {
(0, 1): 10,
(0, 2): 15,
(0, 3): 20,
(1, 2): 35,
(1, 3): 25,
(2, 3): 30
}So, we update it such that the graph is undirected.
# Ensure symmetric distances
for (i, j), dist in list(distances.items()):
if (j, i) not in distances:
distances[j, i] = distSetting Up Variables
We define binary variables for x[i, j] for edges (0 means not connected and 1 means connected).
import dimod
# Binary variables x[i, j]:
# - 1 if the route travels directly from city i to city j
# - 0 otherwise
x = {}
for (i, j) in distances.keys():
x[i, j] = dimod.Binary(f'x[{i},{j}]')Setting Up Objective
We set up the objective using the binary variables and distances.
# Objective: Minimize total distance
objective = 0
for (i, j), dist in distances.items():
objective += dist * x[i, j]Or you can write it more compactly:
# Objective: Minimize total distance
objective = sum(dist * x[i, j] for (i, j), dist in distances.items())Setting Up a CQM Object
We create a CQM object and set the objective to it:
# Create a CQM model
cqm = dimod.ConstrainedQuadraticModel()
cqm.set_objective(objective)Constraint 1: Each city is entered exactly once
For each city j, there should be only one city i that enters into:
# Constraint 1: Each city is entered exactly once
for j in range(num_cities):
# Entries for city j
entries = [x[i, j] for i in range(num_cities) if (i, j) in x and i != j]
if len(entries) > 0:
cqm.add_constraint(sum(entries) == 1, label=f"enter_once_{j}")In this example, we don’t need to check the length of entries but I added there as a protective measure.
Constraint 2: Each city is exited exactly once
For each city i, there should be only one city j that exits into:
# Constraint 2: Each city is exited exactly once
for i in range(num_cities):
# Exits from city i
exits = [x[i, j] for j in range(num_cities) if (i, j) in x and i != j]
if len(exits) > 0:
cqm.add_constraint(sum(exits) == 1, label=f"exit_once_{i}")Constraints 3: Subtour elimination (MTZ constraints)
We set up auxiliary integer variables for ordering cities in the tour:
# Integer variables u[i]: Order in the tour
u = {}
for i in range(num_cities):
u[i] = dimod.Integer(f'u[{i}]', lower_bound=0, upper_bound=num_cities-1)It also ensures that the ordering index is between 0 and the number of cities - 1. We set the starting city as A by the below constraint:
# Explicitly set u[0] = <starting city index>
start_city = cities.index('A')
cqm.add_constraint(u[start_city] == 0, label="StartCityOrder")Then, apply the subtour elimination (MTZ) except for the returning to the start city case:
# Subtour elimination constraints
for i, j in x.keys():
if j != start_city: # return to the starting city
cqm.add_constraint(
u[i] - u[j] + num_cities * x[i, j] <= num_cities - 1,
label=f"MTZ_{i}_{j}"
)Solving with Exact Solver
Since this is a small problem, we’ll use ExactCQMSolver to try all the permutations of variables.
# Exact solver
solver = dimod.ExactCQMSolver()
# Try all combinations of variables
resultset = solver.sample_cqm(cqm)It takes like several seconds depending on how much power you have on your computer. For a TSP problem with n cities:
Binary Decision Variables:
- There are N = n² - n = n(n - 1) binary variables x[i, j], excluding i = j case.
- The total combinations for binary variables is 2ᴺ
Integer Decision Variables:
- There are n integer variables, each with n possible values (0, …, n-1).
- The total combinations for integer variables is nⁿ.
So, ExactCQMSolver will evaluate 2ᴺ ⨉ nⁿ combinations. For n = 4 cities:
- Binary variables: N = 4 ⨉ (4–1) = 12, yeilding 2¹² = 4096 combinations.
- Integer variables: nⁿ = 4⁴ = 256 combinations.
- The total is 4096 ⨉ 256 = 1,048,576 combinations.
Then, the solver will apply the constraints to see if they are feasible or not. So, we can use the is_feasible flag to filter out infeasible results:
# Check if there are any feasible solutions
for sample, energy, is_feasible in resultset.data(['sample',
'energy',
'is_feasible']):
if is_feasible:
# Extract visit order (u[i]) and sort cities by it
route = sorted(cities, key=lambda city: sample[f'u[{cities.index(city)}]'])
print(f"Route: {route} (Total Distance: {energy})")Route: ['A', 'B', 'D', 'C'] (Total Distance: 80.0)
Route: ['A', 'C', 'D', 'B'] (Total Distance: 80.0)
Route: ['A', 'D', 'C', 'B'] (Total Distance: 95.0)
Route: ['A', 'B', 'C', 'D'] (Total Distance: 95.0)
Route: ['A', 'C', 'B', 'D'] (Total Distance: 95.0)
Route: ['A', 'D', 'B', 'C'] (Total Distance: 95.0)We got the same solution as D-Wave NetworkX. The first two are the best feasible solution. The second one is the same route as the second one in the reverse order.
The main takeaway is use ExactCQMSolver only for small problems to ensure your code works as expected. Avoid using it for large problems, as the combinatorial explosion means it may not finish within your lifetime!
We could use LeapHybridCQMSampler instead if you sign up with D-Wave Leap, which gives you a limited access time at no charge. But it may still take a long time, depending on the size and complexity of the problem.
I’ve written other articles that uses Leap solvers (Knapsack, Sudoku, and Max-cut).
QUBO Formulation
Next, let’s try QUBO formulation.
Decision Variables
Define binary variables x[i, j]:
- x[i, j] = 1 if the route includes the edge from city i to city j.
- x[i, j] = 0 otherwise.
- x[i, j] must be 0 when i = j (no self-loop).
Objective Function
The goal is to minimize the total distance traveled:
Contraints as Penalty Terms
Each City is Entered Exactly Once: Every city j must have exactly one incoming edge:
λ₁ is the weight on the penalty.
Each City is Exited Exactly Once: Every city i must have exactly one outgoing edge:
λ₂ is the weight on the penalty.
Sequential Transition Encouragement: A route must transition one city after anoter:
When both x[i, j]=1 and x[j, k]=1, it means the tour includes a transition from city i to j and then from j to k. The product x[i, j] x[j, k] contributes −λ₃ (where λ₃ > 0) to the energy, lowering it (rewarding the configuration).
Since the previous two penalty terms encourage every city is visited exactly once, this term (by rewarding valid transitions) aligns with the overall objective of forming a single, connected tour.
Note: This approach avoids implementing explicit subtour elimination, which would otherwise require additional variables. I didn’t test it with larger or more complex TSPs which may require subtours elimination with additional validation or constraints.
Solution 3: Binary Quadratic Model (BQM)
Here, we’ll use BQM to solve the same TSP.
Setting Up the Symmetric Distances
This part is the same as before.
# Example TSP setup: 4 cities
# cities = [city1, city2, ...]
cities = ['A', 'B', 'C', 'D']
num_cities = len(cities)
# distances = {(city1, city2): distance}
distances = {
(0, 1): 10,
(0, 2): 15,
(0, 3): 20,
(1, 2): 35,
(1, 3): 25,
(2, 3): 30
}
# Ensure symmetric distances
for (i, j), dist in list(distances.items()):
if (j, i) not in distances:
distances[j, i] = distSetting Up a BQM Object and Variables
We define binary variables for x[i, j] for edges (0 means not connected and 1 means connected).
import dimod
# Create the QUBO
qubo = dimod.BinaryQuadraticModel('BINARY')
# Add the objective: Minimize total distance
for (i, j), dist in distances.items():
if i != j:
qubo.add_variable(f'x[{i},{j}]', dist)Constraint 1: Each city is entered exactly once
For each city j, there should be only one city i that enters into:
# Add constraints: Each city is entered exactly once
for j in range(num_cities):
entries = [(f'x[{i},{j}]', 1) for i in range(num_cities) if i != j]
qubo.add_linear_equality_constraint(
entries, constant=-1, lagrange_multiplier=1000,
)langrange_multiplier is λ₁ (the penalty weight).
Constraint 2: Each city is exited exactly once
For each city i, there should be only one city j that exits into:
# Add constraints: Each city is exited exactly once
for i in range(num_cities):
exits = [(f'x[{i},{j}]', 1) for j in range(num_cities) if i != j]
qubo.add_linear_equality_constraint(
exits, constant=-1, lagrange_multiplier=1000,
)langrange_multiplier is λ₂ (the penalty weight).
Constraints 3: Encouraging sequential transitions
This encourages the route to travel from city i to city j and then to city k by adding a reward (negative penalty) for each sequential transition.
# Encourage sequential transitions (i -> j -> k)
for i in range(num_cities):
for j in range(num_cities):
if i != j:
for k in range(num_cities):
if k != i and k != j:
# Encourage the route to travel from i to j and then to k
qubo.add_interaction(
f'x[{i},{j}]',
f'x[{j},{k}]',
-1000
)Since other constraints encourage each city is visited exactly once, this term would encourage forming a single, connected tour.
-1000 is a reward (negative penalty) which encourages sequential transitions.
Solving with Simulated Annealing Sampler
We use SimulatedAnnealingSampler (a classical algorithm) to sample random combinations of binary variables:
# Solve using SimulatedAnnealingSampler
sampler = dimod.SimulatedAnnealingSampler()
sampleset = sampler.sample(qubo, num_reads=100)Let’s look at the first (lowest energy) sample:
sampleset.first.sample(np.float64(-3920.0),
{'x[0,1]': np.int8(0),
'x[0,2]': np.int8(1), # A -> C
'x[0,3]': np.int8(0),
'x[1,0]': np.int8(1), # B -> A
'x[1,2]': np.int8(0),
'x[1,3]': np.int8(0),
'x[2,0]': np.int8(0),
'x[2,1]': np.int8(0),
'x[2,3]': np.int8(1), # C -> D
'x[3,0]': np.int8(0),
'x[3,1]': np.int8(1), # D -> B
'x[3,2]': np.int8(0)})This solution is A → C → D → B, which is correct.
There are 100 samples which can have duplicates, non-optimal solutions, and infeasible ones.
Let’s write a function to decode the results.
def decode_and_normalize(sample, start_city=None):
# Decode the sample into an ordered route and calculate total distance
# Extract the edges from the sample
edges = {i: j for (i, j) in distances.keys() if sample.get(f'x[{i},{j}]') == 1}
if len(edges) != num_cities or len(set(edges.values())) != num_cities:
return None, None # Invalid solution
# Start city index
if start_city is None or start_city not in cities:
start_city = 0
else:
start_city = cities.index(start_city)
# Traverse the route
route, current_city, total_distance = [start_city], start_city, 0
while len(route) < num_cities:
next_city = edges.get(current_city)
if next_city is None or next_city in route:
return None, None # Invalid solution
# Update the total distance and route
total_distance += distances[current_city, next_city]
route.append(next_city)
# Move to the next city
current_city = next_city
# Check if the route completes the cycle
if edges.get(current_city) != start_city:
return None, None # Does not complete the cycle
# Complete the cycle
total_distance += distances[current_city, start_city]
named_route = [cities[i] for i in route]
return named_route, total_distance
# Find unique routes and print the ordered route, total distance, and energy
def print_unique_routes(sampleset, start_city=None):
unique_routes = {}
for sample, energy in sampleset.data(['sample', 'energy']):
# Decode the sample into an ordered route and calculate total distance
named_route, total_distance = decode_and_normalize(sample, start_city=start_city)
if named_route is None:
continue
# Ensure the route is unique
route_key = tuple(named_route)
if route_key not in unique_routes:
unique_routes[route_key] = (total_distance, energy)
# Sort the routes by distance
sorted_routes = sorted(unique_routes.items(), key=lambda item: item[1][0])
# Print sorted unique routes
for route_key, (total_distance, energy) in sorted_routes:
print(f"Ordered Route: {list(route_key)} Total Distance: {total_distance} Energy: {energy}")We can specify the start_city as follows (see start_city = 'B'):
# Print unique routes with starting city 'B'
print_unique_routes(sampleset, start_city='B')The result is as follows (this may change every run):
Ordered Route: ['B', 'A', 'C', 'D'] Total Distance: 80 Energy: -3920.0
Ordered Route: ['B', 'D', 'C', 'A'] Total Distance: 80 Energy: -3920.0
Ordered Route: ['B', 'C', 'D', 'A'] Total Distance: 95 Energy: -3905.0
Ordered Route: ['B', 'C', 'A', 'D'] Total Distance: 95 Energy: -3905.0
Ordered Route: ['B', 'A', 'D', 'C'] Total Distance: 95 Energy: -3905.0
Ordered Route: ['B', 'D', 'A', 'C'] Total Distance: 95 Energy: -3905.0The first two are the best (correct) solution in the reverse order.
Direct QPU
Lastly, we use direct QPU to execute the BQM:
# Use Direct QPU Sampler
from dwave.system import DWaveSampler, EmbeddingComposite
# Use the D-Wave system
sampler = EmbeddingComposite(DWaveSampler())
sampleset = sampler.sample(qubo, num_reads=100)
# Print unique routes with starting city 'B'
print_unique_routes(sampleset, start_city='B')Ordered Route: ['B', 'D', 'C', 'A'] Total Distance: 80 Energy: -3920.0
Ordered Route: ['B', 'A', 'C', 'D'] Total Distance: 80 Energy: -3920.0
Ordered Route: ['B', 'C', 'D', 'A'] Total Distance: 95 Energy: -3905.0
Ordered Route: ['B', 'C', 'A', 'D'] Total Distance: 95 Energy: -3905.0
Ordered Route: ['B', 'A', 'D', 'C'] Total Distance: 95 Energy: -3880.0
Ordered Route: ['B', 'D', 'A', 'C'] Total Distance: 95 Energy: -3875.0
