Quantum Software Development Workshop: Solving TSP with Quantum Computing¶
UPC Postgrado en Ingeniería Cuántica — Quantum Software Development
Dr. José Ramón Martínez Saavedra — March 26, 2026
Overview¶
This notebook walks through a complete Travelling Salesman Problem (TSP) — a canonical combinatorial optimization problem — solved using three different quantum computing frameworks:
- Qiskit 1.0 (IBM) — QAOA with the new Primitives API (Sampler/Estimator)
- PennyLane (Xanadu) — VQE with differentiable quantum computing
- D-Wave Ocean — Simulated annealing (classical analogue of quantum annealing)
Learning objectives¶
- Formulate a combinatorial optimization problem as a QUBO (Quadratic Unconstrained Binary Optimization)
- Map a QUBO to a quantum Hamiltonian
- Solve it using three different quantum paradigms
- Compare results, runtimes, and trade-offs
Prerequisites¶
pip install qiskit qiskit-optimization qiskit-algorithms pennylane dwave-ocean-sdk matplotlib
# Common imports
import numpy as np
import matplotlib.pyplot as plt
from itertools import permutations
import time
np.random.seed(42)
1. Problem Setup: Travelling Salesman Problem (TSP)¶
Given n cities with known pairwise distances, find the shortest route that visits each city exactly once and returns to the starting city.
We use a small instance (3 cities) so that:
- We can verify results by brute force
- The quantum circuits are tractable on a laptop simulator
- The QUBO has $n^2 = 9$ binary variables (= 9 qubits)
# Define cities
n_cities = 3
cities = np.array([
[0.0, 0.0], # City A (depot)
[1.0, 2.0], # City B
[3.0, 1.0], # City C
])
city_labels = ['A', 'B', 'C']
# Compute distance matrix
dist_matrix = np.zeros((n_cities, n_cities))
for i in range(n_cities):
for j in range(n_cities):
dist_matrix[i][j] = np.linalg.norm(cities[i] - cities[j])
print("Distance matrix:")
print(np.round(dist_matrix, 2))
# Visualize the cities
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
ax.scatter(cities[:, 0], cities[:, 1], s=200, c='steelblue', zorder=5)
for i, label in enumerate(city_labels):
ax.annotate(f' {label} ({cities[i,0]:.0f},{cities[i,1]:.0f})',
(cities[i, 0], cities[i, 1]), fontsize=12, fontweight='bold')
# Draw all edges with distances
for i in range(n_cities):
for j in range(i+1, n_cities):
ax.plot([cities[i,0], cities[j,0]], [cities[i,1], cities[j,1]],
'k--', alpha=0.3, linewidth=1)
mid = (cities[i] + cities[j]) / 2
ax.annotate(f'{dist_matrix[i][j]:.2f}', mid, fontsize=9,
ha='center', color='gray')
ax.set_title('TSP Instance: 3 Cities', fontsize=14)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-0.5, 4)
ax.set_ylim(-0.5, 3.5)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Brute-force solution (for verification)
def route_distance(route, dm):
"""Total distance of a closed tour."""
n = len(route)
return sum(dm[route[i]][route[(i+1) % n]] for i in range(n))
best_cost = float('inf')
best_route = None
print("All possible tours (fixing start at A):")
for perm in permutations(range(n_cities)):
cost = route_distance(list(perm), dist_matrix)
labels = ' → '.join(city_labels[p] for p in perm) + f' → {city_labels[perm[0]]}'
marker = ' ← optimal' if cost <= best_cost else ''
if cost < best_cost:
best_cost = cost
best_route = list(perm)
print(f' {labels}: {cost:.4f}{marker}')
print(f'\nOptimal route: {[city_labels[i] for i in best_route]}, cost: {best_cost:.4f}')
2. QUBO Formulation¶
To solve TSP on a quantum computer, we encode it as a QUBO (Quadratic Unconstrained Binary Optimization) problem.
Decision variables¶
$x_{i,t} \in \{0, 1\}$: city $i$ is visited at step $t$
Objective¶
$$\min \sum_{t} \sum_{i,j} d_{ij} \cdot x_{i,t} \cdot x_{j,t+1}$$
Constraints (encoded as penalty terms)¶
- Each city visited exactly once: $\sum_t x_{i,t} = 1 \quad \forall i$
- Each time step has exactly one city: $\sum_i x_{i,t} = 1 \quad \forall t$
These constraints are added to the objective with a penalty multiplier $\lambda$:
$$H = H_{\text{cost}} + \lambda \cdot H_{\text{constraints}}$$
For $n=3$ cities, we need $n^2 = 9$ binary variables (= 9 qubits).
We use qiskit-optimization to build the QUBO automatically:
from qiskit_optimization.applications import Tsp
from qiskit_optimization.converters import QuadraticProgramToQubo
# Build the TSP quadratic program
tsp_app = Tsp(dist_matrix)
qp = tsp_app.to_quadratic_program()
print(qp.prettyprint())
# Convert to QUBO (constraints → penalty terms)
converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)
n_qubits = qubo.get_num_vars()
print(f"QUBO variables (= qubits): {n_qubits}")
print(f"\nQUBO objective:")
print(f" Linear terms: {len(qubo.objective.linear.to_dict())}")
print(f" Quadratic terms: {len(qubo.objective.quadratic.to_dict())}")
print(f" Constant offset: {qubo.objective.constant:.4f}")
3. Approach 1: Qiskit 1.0 — QAOA with Primitives API¶
Qiskit 1.0 (released 2025) introduced a fundamentally new architecture based on Primitives:
- Sampler: returns measurement outcomes (shots)
- Estimator: computes expectation values
The old execute() workflow with Terra/Aer/Ignis/Aqua is deprecated.
QAOA (Quantum Approximate Optimization Algorithm)¶
QAOA alternates between:
- Cost layer: $e^{-i \gamma H_C}$ (encodes the problem)
- Mixer layer: $e^{-i \beta H_M}$ (explores the solution space)
Repeated $p$ times ("reps"), with parameters $(\gamma_1, \beta_1, ..., \gamma_p, \beta_p)$ optimized classically.
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import StatevectorSampler
# --- Step 1: Exact solver (classical baseline, for verification) ---
exact_solver = NumPyMinimumEigensolver()
exact_optimizer = MinimumEigenOptimizer(exact_solver)
result_exact = exact_optimizer.solve(qubo)
route_exact = tsp_app.interpret(result_exact)
dist_exact = route_distance(route_exact, dist_matrix)
print("=== Exact Solver (classical baseline) ===")
print(f"Route: {[city_labels[i] for i in route_exact]}")
print(f"Distance: {dist_exact:.4f}")
print(f"Matches brute force? {abs(dist_exact - best_cost) < 0.01}")
# --- Step 2: QAOA (quantum solver) ---
print("=== QAOA Solver (Qiskit 1.0 Primitives) ===")
print("Running QAOA with StatevectorSampler...")
print("(This may take 1-3 minutes on 9 qubits)\n")
t0 = time.time()
# Key Qiskit 1.0 concept: Primitives replace the old execute() workflow
sampler = StatevectorSampler() # Local statevector-based sampler
qaoa = QAOA(
sampler=sampler,
optimizer=COBYLA(maxiter=200), # Classical optimizer
reps=2, # Number of QAOA layers (p)
)
qaoa_optimizer = MinimumEigenOptimizer(qaoa)
result_qaoa = qaoa_optimizer.solve(qubo)
t_qiskit = time.time() - t0
route_qaoa = tsp_app.interpret(result_qaoa)
valid = all(isinstance(x, (int, np.integer)) for x in route_qaoa)
if valid and len(route_qaoa) == n_cities and len(set(route_qaoa)) == n_cities:
dist_qaoa = route_distance(route_qaoa, dist_matrix)
print(f"Route: {[city_labels[i] for i in route_qaoa]}")
print(f"Distance: {dist_qaoa:.4f}")
print(f"Optimal? {abs(dist_qaoa - best_cost) < 0.01}")
else:
print(f"QAOA returned an invalid tour: {route_qaoa}")
print("This can happen — QAOA is a heuristic. Try increasing reps or maxiter.")
dist_qaoa = None
print(f"Time: {t_qiskit:.1f}s")
print(f"QUBO objective value: {result_qaoa.fval:.4f}")
Key takeaways — Qiskit 1.0¶
StatevectorSampleris the new primitive for local simulation (replacesAer.get_backend('qasm_simulator'))- For real hardware: use
qiskit-ibm-runtimewithSamplerV2connected to IBM Quantum - QAOA is a heuristic — it may not find the optimal solution, especially with few layers
- The classical optimizer (COBYLA) explores the parameter landscape — more iterations help
4. Approach 2: PennyLane — VQE with Differentiable Quantum Computing¶
PennyLane (Xanadu) is the leading framework for differentiable quantum computing — it treats quantum circuits as differentiable functions, enabling gradient-based optimization.
Key differences from Qiskit's QAOA:
- PennyLane computes exact gradients of the quantum circuit (via parameter-shift rule or backpropagation)
- We use a hardware-efficient ansatz (VQE-style) instead of the problem-specific QAOA ansatz
- Native integration with ML frameworks (PyTorch, JAX, TensorFlow)
Mapping QUBO → Ising Hamiltonian¶
PennyLane works with Pauli operators. We map the QUBO to an Ising Hamiltonian:
$$x_i = \frac{1 - Z_i}{2} \quad \Rightarrow \quad H_{\text{QUBO}} \to H_{\text{Ising}}$$
import pennylane as qml
from pennylane import numpy as pnp
# Extract QUBO coefficients
linear = qubo.objective.linear.to_dict()
quadratic = qubo.objective.quadratic.to_dict()
offset = qubo.objective.constant
# Build Ising Hamiltonian from QUBO
# Substitution: x_i = (1 - Z_i) / 2
coeffs_list = []
obs_list = []
const = offset
# Linear terms: v * x_i = v * (1 - Z_i)/2 = v/2 - v/2 * Z_i
for i, v in linear.items():
const += v / 2
coeffs_list.append(-v / 2)
obs_list.append(qml.Z(i))
# Quadratic terms: v * x_i * x_j = v * (1-Z_i)(1-Z_j)/4
# = v/4 * (1 - Z_i - Z_j + Z_i*Z_j)
for (i, j), v in quadratic.items():
if i == j:
const += v / 2
coeffs_list.append(-v / 2)
obs_list.append(qml.Z(i))
else:
const += v / 4
coeffs_list.append(v / 4)
obs_list.append(qml.Z(i) @ qml.Z(j))
coeffs_list.append(-v / 4)
obs_list.append(qml.Z(i))
coeffs_list.append(-v / 4)
obs_list.append(qml.Z(j))
# Include constant as Identity term
coeffs_list.append(const)
obs_list.append(qml.Identity(0))
cost_hamiltonian = qml.Hamiltonian(coeffs_list, obs_list)
cost_hamiltonian = qml.simplify(cost_hamiltonian)
print(f"Ising Hamiltonian: {len(cost_hamiltonian.coeffs)} terms on {n_qubits} qubits")
# VQE with hardware-efficient ansatz
dev = qml.device("default.qubit", wires=n_qubits)
depth = 3 # Ansatz depth
@qml.qnode(dev)
def vqe_cost(params):
"""Hardware-efficient ansatz: layers of RY, RZ rotations + CNOT entanglement."""
for d in range(depth):
# Single-qubit rotations
for i in range(n_qubits):
qml.RY(params[d * n_qubits * 2 + i], wires=i)
qml.RZ(params[d * n_qubits * 2 + n_qubits + i], wires=i)
# Entangling layer (circular)
for i in range(n_qubits - 1):
qml.CNOT(wires=[i, i + 1])
qml.CNOT(wires=[n_qubits - 1, 0])
return qml.expval(cost_hamiltonian)
# Number of variational parameters
n_params = depth * n_qubits * 2
print(f"Ansatz: {depth} layers, {n_params} parameters")
print(f"Circuit depth: ~{depth * 3} gate layers")
# Optimize with Adam (gradient-based — this is PennyLane's strength)
print("=== PennyLane VQE Optimization ===")
print("Using gradient-based Adam optimizer (PennyLane auto-differentiates the circuit)\n")
np.random.seed(42)
params = pnp.array(np.random.uniform(-np.pi, np.pi, n_params), requires_grad=True)
t0 = time.time()
opt = qml.AdamOptimizer(stepsize=0.05)
cost_history = []
best_cost_seen = float('inf')
best_params = params.copy()
n_steps = 200
for step in range(n_steps):
params, cost = opt.step_and_cost(vqe_cost, params)
cost_history.append(float(cost))
if cost < best_cost_seen:
best_cost_seen = float(cost)
best_params = params.copy()
if step % 50 == 0:
print(f" Step {step:3d}: cost = {cost:.4f} (best so far: {best_cost_seen:.4f})")
t_pennylane = time.time() - t0
print(f"\nFinal best cost: {best_cost_seen:.4f}")
print(f"Time: {t_pennylane:.1f}s")
# Plot convergence
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(cost_history, linewidth=1.5, color='darkorange')
ax.axhline(y=best_cost, color='green', linestyle='--', label=f'Optimal = {best_cost:.4f}')
ax.set_xlabel('Optimization step')
ax.set_ylabel('Cost (Hamiltonian expectation)')
ax.set_title('PennyLane VQE Convergence')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Sample the optimized circuit to extract the solution
@qml.qnode(dev)
def sample_circuit(params):
for d in range(depth):
for i in range(n_qubits):
qml.RY(params[d * n_qubits * 2 + i], wires=i)
qml.RZ(params[d * n_qubits * 2 + n_qubits + i], wires=i)
for i in range(n_qubits - 1):
qml.CNOT(wires=[i, i + 1])
qml.CNOT(wires=[n_qubits - 1, 0])
return qml.probs(wires=range(n_qubits))
probs = sample_circuit(best_params)
# Decode top solutions
print("Top 5 measurement outcomes:")
top_indices = np.argsort(probs)[-5:][::-1]
route_pl = None
dist_pl = None
for idx in top_indices:
bs = format(idx, f'0{n_qubits}b')
sol = [int(b) for b in bs]
rm = np.array(sol).reshape(n_cities, n_cities)
route = []
valid = True
for step in range(n_cities):
found = False
for city in range(n_cities):
if rm[city][step] == 1:
route.append(city)
found = True
break
if not found:
valid = False
if valid and len(route) == n_cities and len(set(route)) == n_cities:
dist = route_distance(route, dist_matrix)
labels = [city_labels[i] for i in route]
opt_mark = ' ← optimal' if abs(dist - best_cost) < 0.01 else ''
print(f" |{bs}⟩ p={probs[idx]:.4f} → {labels} dist={dist:.4f}{opt_mark}")
if route_pl is None:
route_pl = route
dist_pl = dist
else:
print(f" |{bs}⟩ p={probs[idx]:.4f} → invalid tour")
Key takeaways — PennyLane¶
- PennyLane auto-differentiates through the quantum circuit — no finite-difference approximation needed
- The
default.qubitdevice simulates the full statevector locally - For real hardware: swap to
qml.device('qiskit.ibmq', ...)orqml.device('amazon.braket', ...) - VQE is more flexible than QAOA (custom ansatz) but has more parameters to optimize
- Gradient-based optimizers (Adam) can be more efficient than gradient-free ones (COBYLA)
5. Approach 3: D-Wave Ocean — Simulated Annealing¶
D-Wave pioneered quantum annealing for optimization. Their Ocean SDK lets you:
- Define a BQM (Binary Quadratic Model) — equivalent to QUBO
- Submit to a quantum annealer or simulated annealer
We use neal.SimulatedAnnealingSampler (runs locally, no D-Wave account needed) to demonstrate the workflow. The same code works with DWaveSampler for real quantum hardware.
Key difference from gate-based approaches¶
- No circuit design needed — you just provide the QUBO/BQM
- No parameter optimization — the annealing schedule is fixed
- Multiple reads — sample many solutions and pick the best one
import dimod
import neal
# Build BQM from the QUBO coefficients
bqm = dimod.BinaryQuadraticModel('BINARY')
# Add linear terms
for i, v in linear.items():
bqm.add_variable(i, v)
# Add quadratic terms (filter self-interactions → they are linear)
for (i, j), v in quadratic.items():
if i != j:
bqm.add_interaction(i, j, v)
else:
bqm.add_variable(i, v)
bqm.offset = offset
print(f"BQM: {bqm.num_variables} variables, {bqm.num_interactions} interactions")
print(f"Offset: {bqm.offset:.4f}")
# Solve with simulated annealing
print("=== D-Wave Ocean: Simulated Annealing ===")
t0 = time.time()
sa_sampler = neal.SimulatedAnnealingSampler()
sampleset = sa_sampler.sample(bqm, num_reads=1000, seed=42)
t_dwave = time.time() - t0
# Best solution
best_sa = sampleset.first
sol = [int(best_sa.sample.get(i, 0)) for i in range(n_qubits)]
rm = np.array(sol).reshape(n_cities, n_cities)
# Decode route
route_dw = []
for step in range(n_cities):
for city in range(n_cities):
if rm[city][step] == 1:
route_dw.append(city)
break
dist_dw = route_distance(route_dw, dist_matrix)
print(f"Route: {[city_labels[i] for i in route_dw]}")
print(f"Distance: {dist_dw:.4f}")
print(f"Energy: {best_sa.energy:.4f}")
print(f"Optimal? {abs(dist_dw - best_cost) < 0.01}")
print(f"Time: {t_dwave:.3f}s")
print(f"Samples taken: {len(sampleset)}")
# Analyze the full sample set
valid_count = 0
optimal_count = 0
for sample in sampleset.samples():
sol = [int(sample.get(i, 0)) for i in range(n_qubits)]
rm = np.array(sol).reshape(n_cities, n_cities)
route = []
for step in range(n_cities):
for city in range(n_cities):
if rm[city][step] == 1:
route.append(city)
break
if len(route) == n_cities and len(set(route)) == n_cities:
valid_count += 1
if abs(route_distance(route, dist_matrix) - best_cost) < 0.01:
optimal_count += 1
print(f"Out of {len(sampleset)} samples:")
print(f" Valid tours: {valid_count} ({100*valid_count/len(sampleset):.1f}%)")
print(f" Optimal tours: {optimal_count} ({100*optimal_count/len(sampleset):.1f}%)")
Key takeaways — D-Wave Ocean¶
- No circuit design: just define the QUBO/BQM and submit
- Extremely fast: simulated annealing finds optimal in milliseconds
- Multiple samples: natural way to explore the solution landscape
- For real hardware: replace
neal.SimulatedAnnealingSampler()withdwave.system.DWaveSampler() - D-Wave Advantage2: 5000+ qubits, but limited to QUBO/Ising problems (no universal computation)
# Summary table
print(f"{'Framework':<20} {'Route':<15} {'Distance':<12} {'Optimal?':<10} {'Time':<10}")
print("=" * 67)
print(f"{'Brute Force':<20} {str([city_labels[i] for i in best_route]):<15} {best_cost:<12.4f} {'YES':<10} {'<0.01s':<10}")
# Qiskit exact
print(f"{'Qiskit (exact)':<20} {str([city_labels[i] for i in route_exact]):<15} {dist_exact:<12.4f} {'YES':<10} {'<1s':<10}")
# Qiskit QAOA
if dist_qaoa is not None:
is_opt = 'YES' if abs(dist_qaoa - best_cost) < 0.01 else 'NO'
print(f"{'Qiskit (QAOA)':<20} {str([city_labels[i] for i in route_qaoa]):<15} {dist_qaoa:<12.4f} {is_opt:<10} {f'{t_qiskit:.1f}s':<10}")
else:
print(f"{'Qiskit (QAOA)':<20} {'invalid':<15} {'N/A':<12} {'NO':<10} {f'{t_qiskit:.1f}s':<10}")
# PennyLane
if dist_pl is not None:
is_opt = 'YES' if abs(dist_pl - best_cost) < 0.01 else 'NO'
print(f"{'PennyLane (VQE)':<20} {str([city_labels[i] for i in route_pl]):<15} {dist_pl:<12.4f} {is_opt:<10} {f'{t_pennylane:.1f}s':<10}")
else:
print(f"{'PennyLane (VQE)':<20} {'no valid':<15} {'N/A':<12} {'NO':<10} {f'{t_pennylane:.1f}s':<10}")
# D-Wave
is_opt = 'YES' if abs(dist_dw - best_cost) < 0.01 else 'NO'
print(f"{'D-Wave (SA)':<20} {str([city_labels[i] for i in route_dw]):<15} {dist_dw:<12.4f} {is_opt:<10} {f'{t_dwave:.3f}s':<10}")
# Visualize the best route found
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
routes = {
'Qiskit QAOA': route_qaoa if (dist_qaoa is not None) else route_exact,
'PennyLane VQE': route_pl if (dist_pl is not None) else best_route,
'D-Wave SA': route_dw,
}
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']
for ax, (name, route), color in zip(axes, routes.items(), colors):
ax.scatter(cities[:, 0], cities[:, 1], s=150, c='steelblue', zorder=5)
for i, label in enumerate(city_labels):
ax.annotate(f' {label}', (cities[i, 0], cities[i, 1]),
fontsize=11, fontweight='bold')
# Draw route
if route is not None:
for k in range(len(route)):
i, j = route[k], route[(k+1) % len(route)]
ax.plot([cities[i,0], cities[j,0]], [cities[i,1], cities[j,1]],
'->', color=color, linewidth=2.5, markersize=8)
dist = route_distance(route, dist_matrix)
ax.set_title(f'{name}\nDist: {dist:.2f}', fontsize=11)
else:
ax.set_title(f'{name}\n(no valid route)', fontsize=11)
ax.set_xlim(-0.5, 4)
ax.set_ylim(-0.5, 3.5)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.suptitle('TSP Solutions: Three Quantum Approaches', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()
7. Discussion & Scaling Considerations¶
What we demonstrated¶
| Aspect | Qiskit QAOA | PennyLane VQE | D-Wave SA |
|---|---|---|---|
| Paradigm | Gate-based (variational) | Gate-based (variational) | Annealing |
| Circuit design | Problem-specific (QAOA) | Generic ansatz (VQE) | No circuit |
| Optimization | Gradient-free (COBYLA) | Gradient-based (Adam) | Physics-based |
| Qubits | n² | n² | n² |
| Strength | IBM hardware ecosystem | Differentiable, flexible | Fast, simple API |
Scaling reality check¶
| Cities | Qubits | Classical brute force | Quantum (current HW) |
|---|---|---|---|
| 3 | 9 | Trivial | Works on simulator |
| 5 | 25 | Easy | Possible on real HW |
| 10 | 100 | Hard | At the edge of current HW |
| 20 | 400 | Very hard | Requires error correction |
| 50 | 2500 | Intractable | Future fault-tolerant QC |
Further reading¶
- arXiv:2502.17725 — VRP on quantum computers (detailed study)
- IBM Quantum Learning — Qiskit tutorials
- PennyLane Demos — QAOA and VQE tutorials
- D-Wave Documentation — Ocean SDK guides
Exercises for the reader¶
- Increase the problem size: Try 4 cities (16 qubits). How does runtime scale for each framework?
- Tune QAOA: Experiment with
reps=1,2,3,4andmaxiter=50,100,200,500. What's the trade-off? - Compare optimizers: Replace COBYLA with SPSA in Qiskit. Replace Adam with GradientDescent in PennyLane.
- Add noise: Use
qiskit.providers.fake_providerto simulate noisy hardware. How do results degrade? - D-Wave vs. classical: Compare D-Wave SA with
scipy.optimize.minimizeon the same QUBO. At what problem size does SA become competitive?
Notebook prepared for the UPC Postgrado en Ingeniería Cuántica, March 2026.
Contact: training@quside.com