Solving Network Problems with D-Wave Ocean: TSP and Graph Coloring

D-Wave quantum annealers are purpose-built for combinatorial optimization. Rather than expressing problems as quantum circuits, you frame them as energy minimization over binary variables, specifically a Quadratic Unconstrained Binary Optimization (QUBO) or a Binary Quadratic Model (BQM). This tutorial covers the physics behind quantum annealing, the mathematical foundations of problem formulation, and several complete worked examples: the Traveling Salesman Problem (TSP), graph coloring, job shop scheduling, and maximum network flow.

Prerequisites

Install the Ocean SDK and set up your D-Wave Leap account:

pip install dwave-ocean-sdk
dwave setup  # authenticate with your Leap API token

You also need NetworkX and NumPy, both of which the Ocean SDK installs as dependencies:

pip install networkx numpy

How Quantum Annealing Works

Before diving into problem formulation, it helps to understand what happens physically inside a D-Wave QPU.

The Hardware

The D-Wave Advantage QPU contains over 5000 superconducting flux qubits arranged in a Pegasus graph topology (specifically P16, with 5627 qubits and roughly 40,000 couplers). Each qubit is a tiny superconducting loop whose persistent current flows either clockwise or counterclockwise, representing the two computational states. Couplers between qubits allow programmable interactions that encode the optimization problem.

The Annealing Process

Quantum annealing solves optimization problems by exploiting the adiabatic theorem of quantum mechanics. The system evolves under a time-dependent Hamiltonian:

$H(s) = (1 - s) \, H_B + s \, H_P$

where $s = s(t)$ is the annealing parameter that increases from 0 to 1 over the annealing time, $H_B$ is the transverse field Hamiltonian (also called the driver Hamiltonian), and $H_P$ is the problem Hamiltonian encoding the optimization objective.

The process follows three stages:

1. Initialization (s = 0): All qubits sit in the ground state of $H_B$, which places every qubit in an equal superposition of its two states. The system “sees” all possible configurations simultaneously through quantum superposition.

2. Annealing (s: 0 to 1): The QPU gradually reduces the transverse field $H_B$ while increasing the problem Hamiltonian $H_P$. If this transition is slow enough relative to the minimum energy gap between ground and first excited states, the adiabatic theorem guarantees the system remains in the instantaneous ground state throughout.

3. Readout (s = 1): The transverse field is fully removed. The system collapses into a classical state that encodes a low-energy (ideally minimum-energy) solution to $H_P$.

A simplified annealing schedule looks like this:

Energy
  |
  |  H_B(s)                H_P(s)
  |  \                      /
  |   \                    /
  |    \                  /
  |     \      gap       /
  |      \    /   \     /
  |       \  /     \   /
  |        \/       \ /
  |        /\        X
  |       /  \      / \
  |      /    \    /   \
  |     /      \  /     \
  |----/--------\/-------\-----> s (0 to 1)
  0                           1

The critical point is where the energy gap between the ground state and first excited state is smallest. If the annealing passes through this point too quickly, the system “jumps” to an excited state and returns a suboptimal solution.

Why Real Annealers Are Not Perfectly Adiabatic

In practice, D-Wave systems use a default annealing time of 20 microseconds, far shorter than what perfect adiabaticity would require for hard problems. Several factors introduce randomness:

• Finite annealing time: The system does not remain perfectly in the ground state, so excited state populations leak in.
• Thermal noise: The QPU operates at approximately 15 millikelvin, but this is not absolute zero. Thermal fluctuations can excite the system, particularly near the minimum gap.
• Control noise: Small imprecisions in programming the qubit biases and coupler strengths add stochastic variation to each anneal.

These effects mean that D-Wave returns a distribution of solutions across multiple reads (called samples or shots). Running 1000 reads and selecting the lowest-energy result is standard practice. The probabilistic nature is a feature, not a bug: it provides a built-in sampling mechanism for exploring the solution landscape.

BQM, QUBO, and Ising Model

The D-Wave Ocean SDK uses three equivalent mathematical representations for optimization problems. Understanding all three and how they relate is essential for effective problem formulation.

QUBO (Quadratic Unconstrained Binary Optimization)

A QUBO minimizes:

$\min_x \; x^T Q \, x$

where $x \in \{0, 1\}^n$ is a vector of binary variables and $Q$ is an $n \times n$ upper-triangular matrix. The diagonal entries $Q_{ii}$ represent linear terms (since $x_i^2 = x_i$ for binary variables), and the off-diagonal entries $Q_{ij}$ represent quadratic interactions between variables $x_i$ and $x_j$.

Ising Model

The Ising model minimizes:

$\min_s \; \sum_{i < j} J_{ij} \, s_i \, s_j + \sum_i h_i \, s_i$

where $s_i \in \{-1, +1\}$ are spin variables, $J_{ij}$ are coupling strengths between spins, and $h_i$ are local magnetic field biases. This is the native representation that D-Wave hardware executes: each qubit physically represents a spin, and couplers between qubits implement the $J_{ij}$ interactions.

BQM (Binary Quadratic Model)

The BQM is Ocean SDK’s general-purpose container that can hold either representation. It stores linear biases, quadratic interactions, an offset constant, and a variable type (BINARY for QUBO variables, SPIN for Ising variables). The BQM is the object you pass to any D-Wave sampler.

Converting Between Representations

The transformation between binary and spin variables is:

$x_i = \frac{s_i + 1}{2}$

This maps $s_i = -1$ to $x_i = 0$ and $s_i = +1$ to $x_i = 1$. The inverse is $s_i = 2x_i - 1$.

Ocean SDK handles all conversions automatically:

import dimod
import numpy as np

# Start from a QUBO matrix
Q = {(0, 0): -1, (1, 1): -1, (0, 1): 2}

# Create BQM from QUBO
bqm = dimod.BinaryQuadraticModel.from_qubo(Q)
print(f"BQM vartype: {bqm.vartype}")
print(f"BQM linear: {dict(bqm.linear)}")
print(f"BQM quadratic: {dict(bqm.quadratic)}")

# Convert to Ising
h, J, offset = bqm.to_ising()
print(f"\nIsing h (biases): {h}")
print(f"Ising J (couplings): {J}")
print(f"Ising offset: {offset}")

# Convert back to QUBO
Q_back, offset_back = bqm.to_qubo()
print(f"\nQUBO (round-trip): {Q_back}")
print(f"QUBO offset: {offset_back}")

# You can also create a BQM directly in spin (Ising) variables
bqm_spin = dimod.BinaryQuadraticModel(
    {"s0": 0.5, "s1": -0.3},          # linear biases h
    {("s0", "s1"): -1.0},              # couplings J
    0.0,                                 # offset
    vartype=dimod.SPIN
)
Q_from_spin, offset = bqm_spin.to_qubo()
print(f"\nSpin BQM converted to QUBO: {Q_from_spin}")

The key takeaway: formulate your problem in whichever representation is most natural. For constraint satisfaction, QUBO variables (0/1) are usually more intuitive. The SDK converts to Ising internally before sending to the QPU.

The Traveling Salesman Problem as a QUBO

The TSP asks: given $N$ cities and pairwise distances, find the shortest route that visits each city exactly once and returns to the start.

Variable Encoding

Use one-hot encoding: define binary variable $x_{i,p}$ where $x_{i,p} = 1$ means city $i$ is visited at position $p$ in the tour. For $N$ cities, this gives $N^2$ binary variables.

For 4 cities (A, B, C, D) at positions (0, 1, 2, 3):

• $x_{A,0} = 1$ means city A is visited first
• $x_{B,2} = 1$ means city B is visited third
• etc.

Constraint Terms

Two hard constraints must be enforced:

1. Each city is visited exactly once: For each city $i$, $\sum_p x_{i,p} = 1$
2. Each position has exactly one city: For each position $p$, $\sum_i x_{i,p} = 1$

Both constraints take the form $(\sum_k y_k - 1)^2 = 0$, which expands to a quadratic penalty:

$P_{\text{city}} = A \sum_i \left(1 - \sum_p x_{i,p}\right)^2$

$P_{\text{position}} = A \sum_p \left(1 - \sum_i x_{i,p}\right)^2$

where $A$ is a penalty coefficient large enough that violating a constraint costs more than any distance reduction.

Objective (Distance Minimization)

The total tour distance is minimized by:

$P_{\text{dist}} = B \sum_{i \neq j} \sum_p d_{ij} \cdot x_{i,p} \cdot x_{j, (p+1) \bmod N}$

where $d_{ij}$ is the distance from city $i$ to city $j$.

TSP Derivation for N=3: Explicit QUBO Matrix

To build intuition, let us work through the complete QUBO construction for $N = 3$ cities. We have 9 binary variables: $x_{0,0}, x_{0,1}, x_{0,2}, x_{1,0}, x_{1,1}, x_{1,2}, x_{2,0}, x_{2,1}, x_{2,2}$.

City constraints (each city visited exactly once):

For city 0: $(x_{0,0} + x_{0,1} + x_{0,2} - 1)^2$

Expanding: $x_{0,0} + x_{0,1} + x_{0,2} - 2x_{0,0}x_{0,1} ...$ Wait, let us be precise. Since $x_k^2 = x_k$ for binary variables:

$(y_0 + y_1 + y_2 - 1)^2 = y_0 + y_1 + y_2 + 2y_0 y_1 + 2y_0 y_2 + 2y_1 y_2 - 2y_0 - 2y_1 - 2y_2 + 1$ $= -y_0 - y_1 - y_2 + 2y_0 y_1 + 2y_0 y_2 + 2y_1 y_2 + 1$

So for each constraint, variables get a linear bias of $-A$ and each pair within the constraint gets a quadratic interaction of $+2A$. There is also a constant offset of $A$ per constraint (6 constraints total = $6A$).

The following code constructs the full 9x9 QUBO matrix for N=3, then verifies that valid tours have the lowest energies:

import numpy as np
from itertools import permutations

# 3-city TSP: distance matrix
N = 3
dist = np.array([
    [0, 5, 8],
    [5, 0, 3],
    [8, 3, 0],
])

A = 100  # penalty coefficient
B = 1    # distance weight

# Build the 9x9 QUBO matrix explicitly
# Variable index: x_{i,p} -> i * N + p
n_vars = N * N
Q = np.zeros((n_vars, n_vars))

def idx(i, p):
    return i * N + p

# City constraints: each city i visited at exactly one position
for i in range(N):
    for p in range(N):
        Q[idx(i, p), idx(i, p)] += -A  # linear term
    for p1 in range(N):
        for p2 in range(p1 + 1, N):
            Q[idx(i, p1), idx(i, p2)] += 2 * A  # cross term

# Position constraints: each position p has exactly one city
for p in range(N):
    for i in range(N):
        Q[idx(i, p), idx(i, p)] += -A  # linear term
    for i1 in range(N):
        for i2 in range(i1 + 1, N):
            Q[idx(i1, p), idx(i2, p)] += 2 * A  # cross term

# Distance objective
for i in range(N):
    for j in range(N):
        if i != j:
            for p in range(N):
                next_p = (p + 1) % N
                r, c = idx(i, p), idx(j, next_p)
                if r <= c:
                    Q[r, c] += B * dist[i][j]
                else:
                    Q[c, r] += B * dist[i][j]

constant_offset = 2 * N * A  # 6 constraints, each contributing A

print("9x9 QUBO matrix Q:")
print(Q)
print(f"\nConstant offset: {constant_offset}")

# Enumerate all 2^9 = 512 binary strings
best_energy = float("inf")
valid_tours = []
all_results = []

for bits in range(2**n_vars):
    x = np.array([(bits >> k) & 1 for k in range(n_vars)], dtype=float)
    energy = x @ Q @ x + constant_offset

    # Check if this is a valid tour
    assignment = x.reshape(N, N)
    row_sums = assignment.sum(axis=1)
    col_sums = assignment.sum(axis=0)
    is_valid = (
        np.all(row_sums == 1) and
        np.all(col_sums == 1)
    )

    all_results.append((energy, bits, is_valid))
    if is_valid:
        # Decode tour
        tour_order = [int(np.argmax(assignment[:, p])) for p in range(N)]
        tour_dist = sum(
            dist[tour_order[p]][tour_order[(p + 1) % N]] for p in range(N)
        )
        valid_tours.append((energy, tour_order, tour_dist))

    if energy < best_energy:
        best_energy = energy

# Sort all results by energy
all_results.sort()

print("\n--- Top 10 lowest-energy configurations ---")
for energy, bits, is_valid in all_results[:10]:
    x = [(bits >> k) & 1 for k in range(n_vars)]
    tag = "VALID" if is_valid else "invalid"
    print(f"  Energy={energy:8.1f}  x={x}  [{tag}]")

print(f"\n--- All {len(valid_tours)} valid tours ---")
valid_tours.sort()
for energy, tour, tour_dist in valid_tours:
    print(f"  Energy={energy:8.1f}  Tour={tour}  Distance={tour_dist}")

# Verify: all valid tours have lower energy than all invalid solutions
max_valid_energy = max(e for e, _, _ in valid_tours)
min_invalid_energy = min(e for e, _, v in all_results if not v)
print(f"\nMax valid tour energy:     {max_valid_energy}")
print(f"Min invalid config energy: {min_invalid_energy}")
print(f"Gap (should be > 0):       {min_invalid_energy - max_valid_energy}")
assert min_invalid_energy > max_valid_energy, "Penalty A is too small!"
print("Verification passed: all valid tours beat all invalid configurations.")

For $N = 3$ with a symmetric distance matrix, there are $3! = 6$ permutations, which correspond to 3 unique tours (each tour appears twice, once forward and once reversed). The output confirms that every valid tour has strictly lower energy than any invalid configuration.

Full QUBO Construction (General N)

import numpy as np
import dimod
from itertools import product

# 4-city example with asymmetric distances
cities = ["A", "B", "C", "D"]
N = len(cities)

# Distance matrix
dist = np.array([
    [0, 2, 9, 10],
    [1, 0, 6, 4],
    [15, 7, 0, 8],
    [6, 3, 12, 0],
])

# Variable naming: x[i][p] -> "x_i_p"
def var(i, p):
    return f"x_{i}_{p}"

# Penalty weight (should be > max tour distance to enforce constraints)
A = 500  # constraint penalty
B = 1    # distance coefficient

Q = {}

def add_to_Q(u, v, coeff):
    if u == v:
        Q[(u, u)] = Q.get((u, u), 0) + coeff
    else:
        key = (min(u, v), max(u, v))
        Q[key] = Q.get(key, 0) + coeff

# Constraint 1: each city visited exactly once
# Expanding (sum_p x[i][p] - 1)^2: linear terms get -A each, cross terms +2A each
for i in range(N):
    for p in range(N):
        add_to_Q(var(i, p), var(i, p), -A)  # linear: -A per variable
    for p1, p2 in product(range(N), repeat=2):
        if p1 < p2:
            add_to_Q(var(i, p1), var(i, p2), 2 * A)  # cross terms

# Constraint 2: each position has exactly one city
for p in range(N):
    for i in range(N):
        add_to_Q(var(i, p), var(i, p), -A)  # linear: -A per variable
    for i1, i2 in product(range(N), repeat=2):
        if i1 < i2:
            add_to_Q(var(i1, p), var(i2, p), 2 * A)

# Objective: minimize total distance
for i in range(N):
    for j in range(N):
        if i != j:
            for p in range(N):
                next_p = (p + 1) % N
                add_to_Q(var(i, p), var(j, next_p), B * dist[i][j])

print(f"QUBO has {len(Q)} non-zero terms for {N} cities ({N**2} variables)")

Solving with LeapHybridSampler

For small instances (N < 8), the D-Wave QPU can solve TSP directly. For larger instances, LeapHybridSampler combines classical and quantum resources:

from dwave.system import LeapHybridSampler

bqm = dimod.BinaryQuadraticModel.from_qubo(Q)

sampler = LeapHybridSampler()
response = sampler.sample(bqm, time_limit=5)  # seconds

best = response.first.sample
energy = response.first.energy

# Decode the solution
tour = [None] * N
valid = True

for i in range(N):
    for p in range(N):
        if best.get(var(i, p), 0) == 1:
            if tour[p] is not None:
                valid = False  # constraint violated
            tour[p] = cities[i]

if valid and None not in tour:
    # Calculate tour distance
    total_dist = 0
    for p in range(N):
        city_i = cities.index(tour[p])
        city_j = cities.index(tour[(p + 1) % N])
        total_dist += dist[city_i][city_j]
    print(f"Tour: {' -> '.join(tour)} -> {tour[0]}")
    print(f"Total distance: {total_dist}")
else:
    print("No valid tour found. Try increasing penalty A or time_limit.")

Penalty Coefficient Selection

Getting penalty weights right is the main challenge in QUBO formulation. If $A$ is too small, the solver ignores constraints and finds low-energy invalid solutions. If $A$ is too large, the solver spends all its energy satisfying constraints and fails to distinguish between feasible solutions of different quality.

Systematic Approach

For a TSP with maximum distance $D_{\max}$ in the distance matrix, a valid tour visits $N$ edges. The worst-case tour cost is $N \cdot D_{\max}$. Setting $A > N \cdot D_{\max}$ guarantees that violating any single constraint costs more than the worst possible tour, so the solver always prefers a valid tour over an invalid configuration with lower distance.

import numpy as np
import dimod

# Distance matrix for 4 cities
dist = np.array([
    [0, 2, 9, 10],
    [1, 0, 6, 4],
    [15, 7, 0, 8],
    [6, 3, 12, 0],
])
N = dist.shape[0]

# Compute the correct penalty coefficient
D_max = dist.max()
A_minimum = N * D_max
A = A_minimum + 1  # add a margin of 1

print(f"Max distance in matrix: {D_max}")
print(f"Worst-case tour cost (upper bound): N * D_max = {N} * {D_max} = {A_minimum}")
print(f"Penalty coefficient A: {A}")

# Show the effect of penalty coefficient choice
def build_tsp_bqm(dist, A, B=1):
    """Build a TSP BQM for a given distance matrix and penalty weight."""
    N = dist.shape[0]
    bqm = dimod.BinaryQuadraticModel("BINARY")

    def var(i, p):
        return f"x_{i}_{p}"

    # City constraints
    for i in range(N):
        for p in range(N):
            bqm.add_variable(var(i, p), -A)
        for p1 in range(N):
            for p2 in range(p1 + 1, N):
                bqm.add_interaction(var(i, p1), var(i, p2), 2 * A)
        bqm.offset += A

    # Position constraints
    for p in range(N):
        for i in range(N):
            bqm.add_variable(var(i, p), -A)
        for i1 in range(N):
            for i2 in range(i1 + 1, N):
                bqm.add_interaction(var(i1, p), var(i2, p), 2 * A)
        bqm.offset += A

    # Distance objective
    for i in range(N):
        for j in range(N):
            if i != j:
                for p in range(N):
                    next_p = (p + 1) % N
                    bqm.add_interaction(var(i, p), var(j, next_p), B * dist[i][j])

    return bqm

# Test with A too small vs. correct A
for test_A in [5, 20, A]:
    bqm = build_tsp_bqm(dist, test_A)
    response = dimod.ExactSolver().sample(bqm)
    best = response.first.sample

    # Check validity
    assignment = np.zeros((N, N))
    for i in range(N):
        for p in range(N):
            assignment[i, p] = best.get(f"x_{i}_{p}", 0)

    row_ok = all(assignment[i, :].sum() == 1 for i in range(N))
    col_ok = all(assignment[:, p].sum() == 1 for p in range(N))
    valid = row_ok and col_ok

    status = "VALID" if valid else "INVALID"
    print(f"\nA={test_A:5.0f}: energy={response.first.energy:10.1f}  [{status}]")
    if valid:
        tour = [int(np.argmax(assignment[:, p])) for p in range(N)]
        tour_dist = sum(dist[tour[p]][tour[(p + 1) % N]] for p in range(N))
        print(f"         Tour: {tour}, Distance: {tour_dist}")
    else:
        print(f"         Assignment row sums: {[int(assignment[i,:].sum()) for i in range(N)]}")

When A is too small (e.g., A=5), you will see that the solver returns “tours” where some cities appear at multiple positions or some positions are empty, because the energy saved by skipping a costly edge outweighs the constraint penalty. With the correct A, all constraints are satisfied and the solver returns the optimal tour.

Graph Coloring as a BQM

Graph coloring asks: assign colors to graph vertices such that no two adjacent vertices share a color, using the minimum number of colors. It appears in scheduling, register allocation, and frequency assignment.

Formulation

Use binary variables $x_{v,c} = 1$ if vertex $v$ is assigned color $c$. With $V$ vertices and $K$ colors, there are $V \times K$ variables.

Constraint 1: each vertex gets exactly one color: $\sum_c x_{v,c} = 1 \quad \forall v$

Constraint 2: adjacent vertices cannot share a color: $x_{u,c} + x_{v,c} \leq 1 \quad \forall (u,v) \in E, \forall c$

Which in penalty form is: $P_{\text{adj}} = A \sum_{(u,v) \in E} \sum_c x_{u,c} \cdot x_{v,c}$

Implementation with dimod

import dimod
import networkx as nx
from itertools import product

# Petersen graph: a well-known graph with chromatic number 3
G = nx.petersen_graph()
K = 3  # number of colors to try

V = list(G.nodes())
colors = list(range(K))

bqm = dimod.BinaryQuadraticModel("BINARY")

A = 10.0  # penalty weight

# Constraint 1: each vertex gets exactly one color
for v in V:
    # (sum_c x[v,c] - 1)^2 penalty
    for c in colors:
        var_name = f"x_{v}_{c}"
        bqm.add_variable(var_name, -A)  # linear: -2A*1 + A (from ^2) = -A
    for c1, c2 in product(colors, repeat=2):
        if c1 < c2:
            bqm.add_interaction(f"x_{v}_{c1}", f"x_{v}_{c2}", 2 * A)
    # Add constant A to offset but dimod handles it automatically
    bqm.offset += A

# Constraint 2: adjacent vertices cannot share the same color
for u, v in G.edges():
    for c in colors:
        bqm.add_interaction(f"x_{u}_{c}", f"x_{v}_{c}", A)

print(f"BQM: {bqm.num_variables} variables, {bqm.num_interactions} interactions")

Solving and Decoding

from dwave.system import LeapHybridSampler
import dimod

# The Petersen graph BQM has 10 vertices * 3 colors = 30 binary variables.
# ExactSolver enumerates 2^30 = ~1 billion states and will hang.
# Use SimulatedAnnealingSampler instead for local testing.
from dwave.samplers import SimulatedAnnealingSampler
sa_sampler = SimulatedAnnealingSampler()
small_response = sa_sampler.sample(bqm, num_reads=200, num_sweeps=2000, seed=42)
best_sample = small_response.first.sample
best_energy = small_response.first.energy

# For production use, the hybrid sampler handles large graphs efficiently:
# from dwave.system import LeapHybridSampler
# sampler = LeapHybridSampler()
# response = sampler.sample(bqm, time_limit=5)
# best_sample = response.first.sample

# Decode coloring
coloring = {}
valid = True
color_names = ["red", "green", "blue", "yellow"]

for v in V:
    assigned = [c for c in colors if best_sample.get(f"x_{v}_{c}", 0) == 1]
    if len(assigned) != 1:
        valid = False
        print(f"Vertex {v}: invalid assignment {assigned}")
    else:
        coloring[v] = assigned[0]

if valid:
    # Check constraint 2
    conflicts = [(u, v) for u, v in G.edges() if coloring[u] == coloring[v]]
    if conflicts:
        print(f"Conflicts found at edges: {conflicts}")
    else:
        print("Valid coloring found!")
        for v, c in coloring.items():
            print(f"  Vertex {v}: {color_names[c]}")
        print(f"Colors used: {len(set(coloring.values()))}")

Finding the Chromatic Number

The chromatic number $\chi(G)$ is the minimum number of colors needed. To find it, try $K = 2, 3, 4, \ldots$ until a valid coloring exists. The Petersen graph has $\chi = 3$: it is not bipartite (so $K = 2$ fails), but a 3-coloring exists.

import dimod
import networkx as nx
from itertools import product

def build_coloring_bqm(G, K, A=10.0):
    """Build a graph coloring BQM for K colors."""
    bqm = dimod.BinaryQuadraticModel("BINARY")
    V = list(G.nodes())
    colors = list(range(K))

    for v in V:
        for c in colors:
            bqm.add_variable(f"x_{v}_{c}", -A)
        for c1, c2 in product(colors, repeat=2):
            if c1 < c2:
                bqm.add_interaction(f"x_{v}_{c1}", f"x_{v}_{c2}", 2 * A)
        bqm.offset += A

    for u, v in G.edges():
        for c in colors:
            bqm.add_interaction(f"x_{u}_{c}", f"x_{v}_{c}", A)

    return bqm

def solve_and_check_coloring(G, K, A=10.0, num_reads=500, num_sweeps=3000, seed=42):
    """Try to K-color graph G. Returns the coloring dict or None.

    Uses SimulatedAnnealingSampler. ExactSolver is impractical for graphs with
    more than ~20 variables (e.g., Petersen with 10 vertices and K=3 colors
    has 30 variables = 2^30 states).
    """
    from dwave.samplers import SimulatedAnnealingSampler
    bqm = build_coloring_bqm(G, K, A)
    response = SimulatedAnnealingSampler().sample(
        bqm, num_reads=num_reads, num_sweeps=num_sweeps, seed=seed
    )
    sample = response.first.sample

    V = list(G.nodes())
    coloring = {}
    for v in V:
        assigned = [c for c in range(K) if sample.get(f"x_{v}_{c}", 0) == 1]
        if len(assigned) != 1:
            return None
        coloring[v] = assigned[0]

    # Check adjacency constraint
    for u, v in G.edges():
        if coloring[u] == coloring[v]:
            return None

    return coloring

# Find the chromatic number of the Petersen graph
G = nx.petersen_graph()
print(f"Petersen graph: {G.number_of_nodes()} vertices, {G.number_of_edges()} edges")
print(f"Is bipartite: {nx.is_bipartite(G)}")

for K in range(2, 6):
    coloring = solve_and_check_coloring(G, K)
    if coloring is not None:
        color_names = ["red", "green", "blue", "yellow", "purple"]
        print(f"\nK={K}: Valid {K}-coloring found! Chromatic number chi(G) = {K}")
        for v in sorted(coloring.keys()):
            print(f"  Vertex {v}: {color_names[coloring[v]]}")

        # Display as adjacency verification
        print("\nEdge verification (all adjacent pairs have different colors):")
        for u, v in list(G.edges())[:5]:
            print(f"  Edge ({u},{v}): {color_names[coloring[u]]} != {color_names[coloring[v]]}")
        print(f"  ... ({G.number_of_edges()} edges total, all verified)")
        break
    else:
        print(f"K={K}: No valid {K}-coloring exists.")

Job Shop Scheduling

Job shop scheduling is a classic combinatorial optimization problem in manufacturing. Given $n$ jobs and $m$ machines, each job consists of a sequence of operations that must be processed on specific machines in a given order, each with a known duration. The goal is to minimize the makespan (the time at which the last operation finishes).

Problem Setup (2 Jobs, 2 Machines)

Consider a minimal example:

• Job 0: Operation 0 on Machine 0 (duration 2), then Operation 1 on Machine 1 (duration 1)
• Job 1: Operation 0 on Machine 1 (duration 2), then Operation 1 on Machine 0 (duration 1)

The makespan lower bound is 3 (the longest single job takes 3 time units), but both jobs share machines, so scheduling conflicts may force a longer makespan.

QUBO Formulation

Define binary variables $x_{j,o,t} = 1$ if operation $o$ of job $j$ starts at time $t$. With a time horizon $T$ (upper bound on makespan), the total number of variables is $n \cdot m \cdot T$.

Constraint 1 (Assignment): Each operation starts at exactly one time: $\sum_t x_{j,o,t} = 1 \quad \forall j, o$

Constraint 2 (Precedence): Within each job, operation $o+1$ cannot start before operation $o$ finishes: $x_{j,o,t_1} + x_{j,o+1,t_2} \leq 1 \quad \text{for all } t_2 < t_1 + d_{j,o}$

where $d_{j,o}$ is the duration of operation $o$ of job $j$.

Constraint 3 (No Overlap): Two operations on the same machine cannot overlap: $x_{j_1,o_1,t_1} + x_{j_2,o_2,t_2} \leq 1 \quad \text{if operations overlap on the same machine}$

import dimod
import numpy as np
from itertools import product

# Problem definition
# jobs[j] = [(machine, duration), ...] in order of operations
jobs = [
    [(0, 2), (1, 1)],  # Job 0: Machine 0 for 2, then Machine 1 for 1
    [(1, 2), (0, 1)],  # Job 1: Machine 1 for 2, then Machine 0 for 1
]

n_jobs = len(jobs)
n_machines = 2
T = 5  # time horizon (upper bound on makespan)

A = 20.0  # constraint penalty

def var(j, o, t):
    return f"x_{j}_{o}_{t}"

bqm = dimod.BinaryQuadraticModel("BINARY")

# Constraint 1: each operation starts at exactly one time
for j in range(n_jobs):
    for o in range(len(jobs[j])):
        times = list(range(T))
        for t in times:
            bqm.add_variable(var(j, o, t), -A)
        for t1, t2 in product(times, repeat=2):
            if t1 < t2:
                bqm.add_interaction(var(j, o, t1), var(j, o, t2), 2 * A)
        bqm.offset += A

# Constraint 2: precedence within each job
for j in range(n_jobs):
    for o in range(len(jobs[j]) - 1):
        machine_o, dur_o = jobs[j][o]
        for t1 in range(T):
            for t2 in range(T):
                # Operation o+1 cannot start before t1 + dur_o
                if t2 < t1 + dur_o:
                    bqm.add_interaction(var(j, o, t1), var(j, o + 1, t2), A)

# Constraint 3: no overlap on same machine
# Build a map: machine -> list of (job, operation, duration)
machine_ops = {}
for j in range(n_jobs):
    for o, (m, d) in enumerate(jobs[j]):
        machine_ops.setdefault(m, []).append((j, o, d))

for m, ops in machine_ops.items():
    for idx1 in range(len(ops)):
        for idx2 in range(idx1 + 1, len(ops)):
            j1, o1, d1 = ops[idx1]
            j2, o2, d2 = ops[idx2]
            for t1 in range(T):
                for t2 in range(T):
                    # Operations overlap if their time windows intersect
                    if not (t2 >= t1 + d1 or t1 >= t2 + d2):
                        bqm.add_interaction(var(j1, o1, t1), var(j2, o2, t2), A)

# Objective: minimize makespan
# Add a small linear cost for later start times of final operations
for j in range(n_jobs):
    last_op = len(jobs[j]) - 1
    _, last_dur = jobs[j][last_op]
    for t in range(T):
        completion_time = t + last_dur
        bqm.add_variable(var(j, last_op, t), 0.1 * completion_time)

print(f"Job Shop BQM: {bqm.num_variables} variables, {bqm.num_interactions} interactions")

# Solve with ExactSolver
response = dimod.ExactSolver().sample(bqm)
sample = response.first.sample

# Decode schedule
print("\nOptimal Schedule:")
makespan = 0
for j in range(n_jobs):
    for o in range(len(jobs[j])):
        machine, duration = jobs[j][o]
        for t in range(T):
            if sample.get(var(j, o, t), 0) == 1:
                end = t + duration
                makespan = max(makespan, end)
                print(f"  Job {j}, Op {o}: Machine {machine}, "
                      f"Start={t}, End={end}")

print(f"\nMakespan: {makespan}")

Maximum Network Flow

Given a directed graph with edge capacities, the maximum flow problem finds the largest total flow from a source node to a sink node while respecting capacity limits and flow conservation at intermediate nodes.

QUBO Formulation

For integer capacities, represent the flow on each edge using binary variables. Define $x_{e,k} = 1$ if the $k$-th unit of flow passes through edge $e$. If edge $e$ has capacity $c_e$, then $k \in \{0, 1, \ldots, c_e - 1\}$, and the total flow on edge $e$ is $f_e = \sum_k x_{e,k}$.

Flow conservation at intermediate nodes (not source or sink): $\sum_{e \in \text{in}(v)} f_e = \sum_{e \in \text{out}(v)} f_e \quad \forall v \notin \{s, t\}$

Since the flow variables are sums of binary variables, the equality constraint becomes: $\left(\sum_{e \in \text{in}(v)} \sum_k x_{e,k} - \sum_{e \in \text{out}(v)} \sum_k x_{e,k}\right)^2 = 0$

Objective: maximize total flow out of source (equivalently, into sink). Since we minimize in QUBO, we negate: $\min -\sum_{e \in \text{out}(s)} \sum_k x_{e,k}$

Capacity constraints are built into the variable encoding: edge $e$ with capacity $c_e$ has exactly $c_e$ binary variables, so the flow cannot exceed $c_e$.

import dimod
import numpy as np

# Small network: 5 nodes, 7 directed edges
# Node 0 = source, Node 4 = sink
edges = [
    (0, 1, 3),  # (from, to, capacity)
    (0, 2, 2),
    (1, 2, 1),
    (1, 3, 2),
    (2, 3, 2),
    (2, 4, 1),
    (3, 4, 3),
]

source = 0
sink = 4
nodes = set()
for u, v, _ in edges:
    nodes.add(u)
    nodes.add(v)
intermediate_nodes = nodes - {source, sink}

A = 20.0  # penalty for flow conservation

bqm = dimod.BinaryQuadraticModel("BINARY")

def flow_var(edge_idx, k):
    return f"f_{edge_idx}_{k}"

# Objective: maximize flow out of source (minimize negative flow)
for idx, (u, v, cap) in enumerate(edges):
    if u == source:
        for k in range(cap):
            bqm.add_variable(flow_var(idx, k), -1.0)  # negative = maximize

# Flow conservation at intermediate nodes
for node in intermediate_nodes:
    # Collect incoming and outgoing flow variable names
    in_vars = []
    out_vars = []
    for idx, (u, v, cap) in enumerate(edges):
        if v == node:
            for k in range(cap):
                in_vars.append(flow_var(idx, k))
        if u == node:
            for k in range(cap):
                out_vars.append(flow_var(idx, k))

    # Penalty: (sum_in - sum_out)^2
    # Let z_i = +1 for in_vars, -1 for out_vars
    all_vars = [(var, +1) for var in in_vars] + [(var, -1) for var in out_vars]

    for i, (vi, si) in enumerate(all_vars):
        # Linear term: A * si * si * x_i = A * x_i (since si^2 = 1, x_i^2 = x_i)
        bqm.add_variable(vi, A * si * si)
        # But for (sum z_i x_i)^2, the linear part from x_i^2 is A * 1
        # and we need to subtract 2 from cross with itself... let us just expand properly
        pass

    # Clear and redo: expand (sum_i s_i x_i)^2 = sum_i x_i + 2 * sum_{i<j} s_i s_j x_i x_j
    # since s_i^2 = 1 and x_i^2 = x_i for binary
    # Reset the linear terms we just added
    for vi, si in all_vars:
        bqm.add_variable(vi, -A * si * si)  # undo the above

    for i, (vi, si) in enumerate(all_vars):
        bqm.add_variable(vi, A)  # linear from x_i^2 = x_i
    for i in range(len(all_vars)):
        for j in range(i + 1, len(all_vars)):
            vi, si = all_vars[i]
            vj, sj = all_vars[j]
            bqm.add_interaction(vi, vj, 2 * A * si * sj)

print(f"Max Flow BQM: {bqm.num_variables} variables, {bqm.num_interactions} interactions")

# Solve
response = dimod.ExactSolver().sample(bqm)
sample = response.first.sample

# Decode flows
print("\nEdge flows:")
total_out_source = 0
for idx, (u, v, cap) in enumerate(edges):
    flow = sum(sample.get(flow_var(idx, k), 0) for k in range(cap))
    print(f"  Edge ({u} -> {v}): flow={flow}/{cap}")
    if u == source:
        total_out_source += flow

print(f"\nMaximum flow: {total_out_source}")

# Verify flow conservation
for node in intermediate_nodes:
    flow_in = 0
    flow_out = 0
    for idx, (u, v, cap) in enumerate(edges):
        flow = sum(sample.get(flow_var(idx, k), 0) for k in range(cap))
        if v == node:
            flow_in += flow
        if u == node:
            flow_out += flow
    status = "OK" if flow_in == flow_out else "VIOLATED"
    print(f"  Node {node}: in={flow_in}, out={flow_out} [{status}]")

D-Wave Pegasus Topology and Embedding

When you submit a problem to the D-Wave QPU, the logical QUBO variables must be mapped onto physical qubits. This process, called minor embedding, is necessary because the Pegasus hardware graph is sparse: each qubit connects to at most 15 others, while your QUBO may have interactions between any pair of variables.

How Embedding Works

When two logical variables interact in the QUBO but their corresponding physical qubits are not directly connected in Pegasus, the embedder creates a chain: a group of connected physical qubits that all represent the same logical variable. Strong ferromagnetic couplings within the chain force all qubits in the chain to agree, effectively acting as a single logical variable.

The minorminer library, included with the Ocean SDK, handles this automatically. Here is how to inspect the embedding for a small problem:

import dimod
import networkx as nx
import minorminer
from dwave.system import DWaveSampler

# Build a small TSP QUBO (3 cities) to embed
N = 3
dist = [[0, 5, 8], [5, 0, 3], [8, 3, 0]]
A = 100

bqm = dimod.BinaryQuadraticModel("BINARY")
for i in range(N):
    for p in range(N):
        bqm.add_variable(f"x_{i}_{p}", -A)
    for p1 in range(N):
        for p2 in range(p1 + 1, N):
            bqm.add_interaction(f"x_{i}_{p1}", f"x_{i}_{p2}", 2 * A)
    bqm.offset += A
for p in range(N):
    for i in range(N):
        bqm.add_variable(f"x_{i}_{p}", -A)
    for i1 in range(N):
        for i2 in range(i1 + 1, N):
            bqm.add_interaction(f"x_{i1}_{p}", f"x_{i2}_{p}", 2 * A)
    bqm.offset += A
for i in range(N):
    for j in range(N):
        if i != j:
            for p in range(N):
                next_p = (p + 1) % N
                bqm.add_interaction(f"x_{i}_{p}", f"x_{j}_{next_p}", dist[i][j])

# Get the logical problem graph
source_graph = dimod.to_networkx_graph(bqm)

# Use the Pegasus topology (P16 as used on Advantage)
# For offline work, generate the ideal Pegasus graph:
import dwave_networkx as dnx
pegasus = dnx.pegasus_graph(16)

# Find an embedding
embedding = minorminer.find_embedding(source_graph, pegasus)

if embedding:
    print(f"Embedding found for {len(embedding)} logical variables")
    chain_lengths = {var: len(chain) for var, chain in embedding.items()}
    print(f"Chain lengths: {chain_lengths}")
    print(f"Max chain length: {max(chain_lengths.values())}")
    print(f"Total physical qubits used: {sum(chain_lengths.values())}")
else:
    print("No embedding found. Problem may be too large for this topology.")

Why Chain Length Matters

Long chains introduce two problems:

1. Chain breaks: If any qubit in a chain disagrees with the rest, the chain “breaks” and the logical variable has no clear value. The sampler must choose how to resolve this (majority vote is the default). Longer chains break more often.

2. Reduced precision: The strong ferromagnetic couplings needed to hold chains together consume part of the QPU’s limited dynamic range, leaving less room for the problem’s own biases and couplings.

A rule of thumb: if your maximum chain length exceeds 7 or 8, solution quality degrades noticeably. For problems requiring very long chains, the LeapHybridSampler is a better choice because it handles embedding and decomposition internally.

Solver Comparison: Simulated Annealing, QPU, and Hybrid

The Ocean SDK provides three main solver types, each suited to different scenarios.

SimulatedAnnealingSampler (Classical)

This is a classical heuristic that mimics the annealing process in software. It requires no D-Wave account, runs locally, and is ideal for development and testing.

import dimod
import networkx as nx
from itertools import product
from dwave.samplers import SimulatedAnnealingSampler

# Build graph coloring BQM for a 5-cycle with 3 colors
G = nx.cycle_graph(5)
K = 3
A = 10.0

bqm = dimod.BinaryQuadraticModel("BINARY")
for v in G.nodes():
    for c in range(K):
        bqm.add_variable(f"x_{v}_{c}", -A)
    for c1, c2 in product(range(K), repeat=2):
        if c1 < c2:
            bqm.add_interaction(f"x_{v}_{c1}", f"x_{v}_{c2}", 2 * A)
    bqm.offset += A
for u, v in G.edges():
    for c in range(K):
        bqm.add_interaction(f"x_{u}_{c}", f"x_{v}_{c}", A)

# SimulatedAnnealingSampler with tunable parameters
sa_sampler = SimulatedAnnealingSampler()
response = sa_sampler.sample(
    bqm,
    num_reads=100,       # number of independent annealing runs
    num_sweeps=1000,     # sweeps per run (more = better but slower)
    beta_range=[0.1, 3], # inverse temperature range [start, end]
    seed=42,             # for reproducibility
)

print(f"Best energy: {response.first.energy}")
print(f"Number of samples: {len(response)}")

# Count how many samples found valid colorings
valid_count = 0
for sample, energy in zip(response.samples(), response.data_vectors["energy"]):
    coloring = {}
    ok = True
    for v in G.nodes():
        assigned = [c for c in range(K) if sample.get(f"x_{v}_{c}", 0) == 1]
        if len(assigned) != 1:
            ok = False
            break
        coloring[v] = assigned[0]
    if ok:
        conflicts = [1 for u, v in G.edges() if coloring[u] == coloring[v]]
        if not conflicts:
            valid_count += 1

print(f"Valid colorings found: {valid_count}/{len(response)}")

DWaveSampler (Direct QPU)

This sends the problem directly to D-Wave hardware. Best for small, dense problems that embed into Pegasus without long chains:

from dwave.system import DWaveSampler, EmbeddingComposite

# EmbeddingComposite automatically handles the minor embedding
qpu_sampler = EmbeddingComposite(DWaveSampler())

response = qpu_sampler.sample(
    bqm,
    num_reads=1000,           # number of annealing shots
    annealing_time=20,        # microseconds (default 20, max 2000)
    chain_strength=None,      # auto-scale (or set manually)
)

# Check for chain breaks
for datum in response.data():
    if datum.chain_break_fraction > 0:
        print(f"Chain break fraction: {datum.chain_break_fraction:.3f}")
        break

print(f"Best energy: {response.first.energy}")
print(f"Timing: {response.info['timing']}")

LeapHybridSampler (Hybrid Classical-Quantum)

This is the most versatile option. It decomposes large problems, applies classical preprocessing, and uses the QPU for the hardest subproblems:

from dwave.system import LeapHybridSampler

hybrid_sampler = LeapHybridSampler()

response = hybrid_sampler.sample(
    bqm,
    time_limit=5,  # minimum 3 seconds, longer = potentially better solutions
)

print(f"Best energy: {response.first.energy}")

# LeapHybridSampler returns a single best sample by default
print(f"Number of samples: {len(response)}")

When to Use Which

Scenario	Recommended Solver
Development and testing	SimulatedAnnealingSampler
Small problems (< 150 logical variables)	DWaveSampler via EmbeddingComposite
Large problems (hundreds to millions of variables)	LeapHybridSampler
Benchmarking against classical	SimulatedAnnealingSampler + LeapHybridSampler
Problems with dense connectivity	LeapHybridSampler (avoids long chains)

The LeapHybridSampler time_limit parameter controls the total runtime. The minimum is 3 seconds for most problems. For very large BQMs (over 100,000 variables), increasing time_limit to 10 or 20 seconds often improves solution quality.

Post-Processing and Solution Polishing

D-Wave solutions are often near-optimal but not exactly optimal, especially for problems that require long chains or have complex constraint structures. Classical post-processing can close the gap efficiently.

Steepest Descent (Greedy Local Search)

The SteepestDescentSolver from dwave.preprocessing takes a D-Wave sample and iteratively flips variables that reduce energy. It runs in polynomial time and often finds the nearest local minimum:

import dimod
import numpy as np
from dwave.samplers import SimulatedAnnealingSampler, SteepestDescentSolver

# Build a TSP BQM (reusing the 4-city example)
cities = ["A", "B", "C", "D"]
N = len(cities)
dist = np.array([
    [0, 2, 9, 10],
    [1, 0, 6, 4],
    [15, 7, 0, 8],
    [6, 3, 12, 0],
])

def var(i, p):
    return f"x_{i}_{p}"

A = 61  # N * D_max + 1 = 4 * 15 + 1
bqm = dimod.BinaryQuadraticModel("BINARY")

for i in range(N):
    for p in range(N):
        bqm.add_variable(var(i, p), -A)
    for p1 in range(N):
        for p2 in range(p1 + 1, N):
            bqm.add_interaction(var(i, p1), var(i, p2), 2 * A)
    bqm.offset += A

for p in range(N):
    for i in range(N):
        bqm.add_variable(var(i, p), -A)
    for i1 in range(N):
        for i2 in range(i1 + 1, N):
            bqm.add_interaction(var(i1, p), var(i2, p), 2 * A)
    bqm.offset += A

for i in range(N):
    for j in range(N):
        if i != j:
            for p in range(N):
                next_p = (p + 1) % N
                bqm.add_interaction(var(i, p), var(j, next_p), dist[i][j])

# Step 1: Get initial samples from simulated annealing
sa = SimulatedAnnealingSampler()
raw_response = sa.sample(bqm, num_reads=50, num_sweeps=500, seed=42)

print(f"Raw SA best energy: {raw_response.first.energy}")

# Step 2: Polish with steepest descent
sd = SteepestDescentSolver()
polished_response = sd.sample(bqm, initial_states=raw_response)

print(f"Polished best energy: {polished_response.first.energy}")
print(f"Energy improvement: {raw_response.first.energy - polished_response.first.energy}")

Tabu Search

For more thorough post-processing, TabuSampler implements tabu search, which avoids revisiting recently explored states. It explores more of the solution space than steepest descent at the cost of longer runtime:

from dwave.samplers import TabuSampler

tabu = TabuSampler()
tabu_response = tabu.sample(
    bqm,
    initial_states=raw_response,
    tenure=5,        # how many recent moves to forbid
    num_reads=10,    # independent tabu searches
    timeout=1000,    # milliseconds per read
)

print(f"Tabu best energy: {tabu_response.first.energy}")

# Compare all three
print("\nComparison:")
print(f"  Raw SA:          {raw_response.first.energy:.1f}")
print(f"  Steepest Descent: {polished_response.first.energy:.1f}")
print(f"  Tabu Search:      {tabu_response.first.energy:.1f}")

A practical workflow combines all three: use SimulatedAnnealingSampler (or DWaveSampler) for the initial solution, then apply SteepestDescentSolver for fast local improvement, and finally use TabuSampler for deeper refinement if time permits.

Cost Comparison and Benchmarking

How does D-Wave’s hybrid solver compare to purely classical approaches? The answer depends heavily on problem size and structure.

Benchmarking TSP at Different Scales

import time
import numpy as np
import dimod
from itertools import permutations
from dwave.samplers import SimulatedAnnealingSampler

def random_distance_matrix(n, seed=42):
    """Generate a random symmetric distance matrix."""
    rng = np.random.RandomState(seed)
    d = rng.randint(1, 20, size=(n, n))
    d = (d + d.T) // 2
    np.fill_diagonal(d, 0)
    return d

def build_tsp_bqm(dist):
    """Build TSP BQM with automatic penalty scaling."""
    N = dist.shape[0]
    D_max = dist.max()
    A = N * D_max + 1

    bqm = dimod.BinaryQuadraticModel("BINARY")

    def var(i, p):
        return f"x_{i}_{p}"

    for i in range(N):
        for p in range(N):
            bqm.add_variable(var(i, p), -A)
        for p1 in range(N):
            for p2 in range(p1 + 1, N):
                bqm.add_interaction(var(i, p1), var(i, p2), 2 * A)
        bqm.offset += A

    for p in range(N):
        for i in range(N):
            bqm.add_variable(var(i, p), -A)
        for i1 in range(N):
            for i2 in range(i1 + 1, N):
                bqm.add_interaction(var(i1, p), var(i2, p), 2 * A)
        bqm.offset += A

    for i in range(N):
        for j in range(N):
            if i != j:
                for p in range(N):
                    next_p = (p + 1) % N
                    bqm.add_interaction(var(i, p), var(j, next_p), dist[i][j])

    return bqm

def brute_force_tsp(dist):
    """Find optimal tour by exhaustive enumeration."""
    N = dist.shape[0]
    best_dist = float("inf")
    best_tour = None
    # Fix city 0 at position 0 to reduce permutations
    for perm in permutations(range(1, N)):
        tour = [0] + list(perm)
        d = sum(dist[tour[i]][tour[(i + 1) % N]] for i in range(N))
        if d < best_dist:
            best_dist = d
            best_tour = tour
    return best_tour, best_dist

def decode_tsp_sample(sample, N):
    """Decode a TSP sample into a tour and distance."""
    assignment = np.zeros((N, N))
    for i in range(N):
        for p in range(N):
            assignment[i, p] = sample.get(f"x_{i}_{p}", 0)

    row_ok = all(assignment[i, :].sum() == 1 for i in range(N))
    col_ok = all(assignment[:, p].sum() == 1 for p in range(N))
    if not (row_ok and col_ok):
        return None, None
    tour = [int(np.argmax(assignment[:, p])) for p in range(N)]
    return tour, None  # distance computed separately

# Benchmark for N = 5, 8, 10
for N in [5, 8, 10]:
    dist = random_distance_matrix(N)
    print(f"\n{'='*50}")
    print(f"TSP with N={N} cities ({N**2} binary variables)")
    print(f"{'='*50}")

    # Brute force (only feasible for small N)
    if N <= 10:
        t0 = time.time()
        opt_tour, opt_dist = brute_force_tsp(dist)
        bf_time = time.time() - t0
        print(f"\nBrute force: distance={opt_dist}, time={bf_time:.3f}s")
    else:
        opt_dist = None
        print("\nBrute force: skipped (too large)")

    # Simulated annealing
    bqm = build_tsp_bqm(dist)
    sa = SimulatedAnnealingSampler()
    t0 = time.time()
    sa_response = sa.sample(bqm, num_reads=200, num_sweeps=2000, seed=42)
    sa_time = time.time() - t0

    sa_best = sa_response.first.sample
    sa_tour, _ = decode_tsp_sample(sa_best, N)
    if sa_tour:
        sa_dist = sum(dist[sa_tour[i]][sa_tour[(i + 1) % N]] for i in range(N))
        ratio = sa_dist / opt_dist if opt_dist else float("nan")
        print(f"SA: distance={sa_dist}, ratio={ratio:.3f}, time={sa_time:.3f}s")
    else:
        print(f"SA: no valid tour found, time={sa_time:.3f}s")

    # LeapHybridSampler (uncomment when you have a Leap account)
    # from dwave.system import LeapHybridSampler
    # hybrid = LeapHybridSampler()
    # t0 = time.time()
    # hybrid_response = hybrid.sample(bqm, time_limit=5)
    # hybrid_time = time.time() - t0
    # hybrid_best = hybrid_response.first.sample
    # hybrid_tour, _ = decode_tsp_sample(hybrid_best, N)
    # if hybrid_tour:
    #     hybrid_dist = sum(
    #         dist[hybrid_tour[i]][hybrid_tour[(i+1) % N]] for i in range(N)
    #     )
    #     ratio = hybrid_dist / opt_dist if opt_dist else float("nan")
    #     print(f"Hybrid: distance={hybrid_dist}, ratio={ratio:.3f}, "
    #           f"time={hybrid_time:.3f}s")

Honest Assessment

For TSP at small scales (N < 15), classical solvers (including branch-and-bound, Concorde, and even brute force) outperform D-Wave both in solution quality and runtime. The QUBO formulation overhead is substantial: N cities require N^2 binary variables, and embedding those into Pegasus adds even more physical qubits.

Where D-Wave’s hybrid solver shows promise is on large, structured optimization problems (thousands to millions of variables) where exact classical solvers are impractical and heuristics struggle to escape local minima. The quantum component provides a different exploration mechanism that sometimes finds solutions classical heuristics miss.

For problems below roughly 200 variables, SimulatedAnnealingSampler is typically sufficient. For problems between 200 and 100,000 variables, LeapHybridSampler offers a convenient, competitive option. Above that, problem decomposition strategies (which Leap handles internally) become essential.

Common Mistakes

Working with D-Wave and QUBO formulations involves several subtle pitfalls. Here are the most common mistakes and how to avoid them.

1. Penalty Coefficient Too Small

If the constraint penalty $A$ is smaller than the maximum possible objective value, the solver finds “solutions” that violate constraints because doing so reduces total energy. Always compute the penalty systematically:

# Wrong: arbitrary small penalty
A = 10  # too small if distances can sum to more than 10

# Right: derive from problem parameters
D_max = dist.max()
A = N * D_max + 1

2. Using QPU Directly for Large Problems

The D-Wave Advantage QPU has roughly 5600 qubits, but embedding overhead typically reduces the effective capacity to 100-200 logical variables. Submitting a 500-variable QUBO to DWaveSampler causes an embedding failure or extremely long chains:

# Wrong: sending a large problem to the QPU
from dwave.system import DWaveSampler, EmbeddingComposite
# This will likely fail or produce poor results for large BQMs
# response = EmbeddingComposite(DWaveSampler()).sample(large_bqm)

# Right: use LeapHybridSampler for large problems
from dwave.system import LeapHybridSampler
response = LeapHybridSampler().sample(large_bqm, time_limit=10)

3. Incorrect Variable Decoding

A common bug is constructing variable names differently during QUBO building and solution decoding. Always use the same naming function:

# Define once and reuse
def var(i, p):
    return f"x_{i}_{p}"

# Building phase
bqm.add_variable(var(i, p), bias)

# Decoding phase (use the SAME function)
for i in range(N):
    for p in range(N):
        if sample[var(i, p)] == 1:  # consistent naming
            tour[p] = i

4. Expecting Deterministic Results

Quantum annealing is inherently probabilistic. Running the same problem twice produces different samples. Always take multiple reads and select the best:

# Wrong: single read, assume the result is optimal
response = sampler.sample(bqm, num_reads=1)

# Right: multiple reads, pick the best
response = sampler.sample(bqm, num_reads=1000)
best = response.first  # lowest energy across all reads
print(f"Best energy: {best.energy}")
print(f"Number of distinct solutions: {len(response)}")

5. Ignoring Chain Breaks

When using DWaveSampler, always check for chain breaks. A high chain break fraction indicates that the embedding is strained and solutions may be unreliable:

# After sampling with DWaveSampler
for datum in response.data():
    if datum.chain_break_fraction > 0.1:
        print(f"WARNING: {datum.chain_break_fraction:.1%} chain break fraction")
        print("Consider: increasing chain_strength, using LeapHybridSampler,")
        print("or reducing problem size.")

6. Using ExactSolver for Large Problems

The ExactSolver enumerates all $2^n$ states, which is useful for verification but impractical beyond about 20 variables. For 30 variables, it evaluates over one billion states:

# This is fine (9 variables = 512 states)
response = dimod.ExactSolver().sample(small_bqm)

# This will hang or crash (40 variables = 1 trillion states)
# response = dimod.ExactSolver().sample(large_bqm)  # DO NOT DO THIS

# For medium problems, use SimulatedAnnealingSampler instead
from dwave.samplers import SimulatedAnnealingSampler
response = SimulatedAnnealingSampler().sample(large_bqm, num_reads=100)

Next Steps

This tutorial covers the core workflow for solving network optimization problems with D-Wave Ocean. To go further:

• Constrained Quadratic Models (CQMs): Ocean’s ConstrainedQuadraticModel class lets you express constraints directly without manual penalty tuning. The LeapHybridCQMSampler handles the conversion internally.
• Discrete Quadratic Models (DQMs): For problems with discrete (non-binary) variables, DQMs avoid the overhead of one-hot encoding.
• Problem decomposition: For very large problems, explore dwave.preprocessing tools like FixVariablesComposite (which fixes variables that can be determined classically) and CutOffsComposite (which prunes weak interactions).
• Real QPU experiments: If you have a D-Wave Leap account, try running the examples on actual hardware and compare the solution distributions to the classical simulators.