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:

  1. Qiskit 1.0 (IBM) — QAOA with the new Primitives API (Sampler/Estimator)
  2. PennyLane (Xanadu) — VQE with differentiable quantum computing
  3. 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
In [ ]:
# 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)
In [ ]:
# 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))
In [ ]:
# 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()
In [ ]:
# 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)¶

  1. Each city visited exactly once: $\sum_t x_{i,t} = 1 \quad \forall i$
  2. 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:

In [ ]:
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())
In [ ]:
# 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:

  1. Cost layer: $e^{-i \gamma H_C}$ (encodes the problem)
  2. 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.

In [ ]:
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}")
In [ ]:
# --- 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¶

  • StatevectorSampler is the new primitive for local simulation (replaces Aer.get_backend('qasm_simulator'))
  • For real hardware: use qiskit-ibm-runtime with SamplerV2 connected 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}}$$

In [ ]:
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")
In [ ]:
# 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")
In [ ]:
# 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")
In [ ]:
# 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()
In [ ]:
# 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.qubit device simulates the full statevector locally
  • For real hardware: swap to qml.device('qiskit.ibmq', ...) or qml.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:

  1. Define a BQM (Binary Quadratic Model) — equivalent to QUBO
  2. 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
In [ ]:
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}")
In [ ]:
# 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)}")
In [ ]:
# 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() with dwave.system.DWaveSampler()
  • D-Wave Advantage2: 5000+ qubits, but limited to QUBO/Ising problems (no universal computation)

6. Comparison of Results¶

Let's compare the three approaches side by side.

In [ ]:
# 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}")
In [ ]:
# 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¶

  1. Increase the problem size: Try 4 cities (16 qubits). How does runtime scale for each framework?
  2. Tune QAOA: Experiment with reps=1,2,3,4 and maxiter=50,100,200,500. What's the trade-off?
  3. Compare optimizers: Replace COBYLA with SPSA in Qiskit. Replace Adam with GradientDescent in PennyLane.
  4. Add noise: Use qiskit.providers.fake_provider to simulate noisy hardware. How do results degrade?
  5. D-Wave vs. classical: Compare D-Wave SA with scipy.optimize.minimize on 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