D-Wave / Volkswagen: Quantum Annealing for Factory Job Shop Scheduling

The Manufacturing Scheduling Problem

Modern automotive body shops operate dozens to hundreds of robotic welding stations. Each car body (a job) must pass through a defined sequence of stations (operations), and each station can handle only one body at a time. The question is: in what order should bodies be processed at each station to minimize overall completion time (makespan)?

This is the Job Shop Scheduling Problem (JSSP), one of the canonical NP-hard combinatorial optimization problems. For a factory with n jobs and m machines, the number of possible schedules is (n!)^m. A plant with 20 car bodies and 30 welding stations has roughly 10^84 candidate schedules.

Volkswagen currently solves this with a combination of constraint programming (IBM CPLEX) for small instances and hand-tuned heuristics for production scheduling. Both approaches struggle with the dynamic replanning needed when a robot goes offline or a new order arrives. Quantum annealing was evaluated as a fast approximate solver that could respond to disruptions in under one second.

Quantum Annealing vs Gate-Based Quantum Computing

D-Wave’s approach is fundamentally different from IBM, IonQ, or Google hardware. Rather than applying sequences of logic gates, D-Wave implements quantum annealing: a physical process that searches for the lowest-energy state of a programmable Ising Hamiltonian.

The processor starts with all qubits in a superposition (quantum tunneling is active), then slowly reduces the tunneling field while increasing the strength of programmable qubit couplings. If this is done slowly enough, the system should settle into the ground state of the target Hamiltonian, which encodes the solution to an optimization problem.

Key differences from gate-based computing:

Property	Gate-based (IBM, IonQ)	Quantum annealing (D-Wave)
Qubit count (2024)	127-1121 qubits	5000+ qubits
Programmability	Universal: any quantum algorithm	Specialized: QUBO/Ising problems only
Coherence time	Microseconds to milliseconds	Microseconds (annealing schedule)
Error correction	Possible in principle	Not applicable
Practical use today	NISQ experiments, benchmarking	Combinatorial optimization heuristic

D-Wave does not claim quantum speedup in the theoretical sense. It positions the Advantage processor as a highly parallel heuristic that often finds good solutions faster than simulated annealing or classical branch-and-bound within a fixed time budget.

QUBO Formulation of Job Shop Scheduling

To run on D-Wave, the scheduling problem must be expressed as a Quadratic Unconstrained Binary Optimization (QUBO) problem. Each binary variable x_{i,k,t} = 1 means job i is processed on machine k starting at time slot t.

The objective and constraints become energy terms:

H = H_objective + lambda_1 * H_one_slot + lambda_2 * H_precedence + lambda_3 * H_capacity

• H_objective: minimize makespan (latest completion time)
• H_one_slot: each job-machine pair occupies exactly one time slot
• H_precedence: operations for each job follow the required sequence
• H_capacity: each machine handles at most one job at a time

The lambda values are penalty weights that must be tuned so constraints are satisfied while the objective is minimized.

Ocean SDK Implementation

D-Wave’s Ocean SDK handles the translation from a high-level problem description to the Ising Hamiltonian required by the hardware.

import dimod
import dwave.inspector
from dwave.system import DWaveSampler, EmbeddingComposite
from collections import defaultdict
import itertools

# ---- Define the problem ----
# jobs[j] = list of (machine, processing_time) pairs in order
jobs = {
    0: [(0, 3), (1, 2), (2, 2)],
    1: [(0, 2), (2, 1), (1, 4)],
    2: [(1, 4), (0, 1), (2, 3)],
}

n_jobs = len(jobs)
n_machines = 3
max_time = 12  # discretized time horizon

# Binary variables: x[j][m][t] = 1 if job j starts on machine m at time t
# Flatten into QUBO variable indices
def var(j, m, t):
    return f"x_{j}_{m}_{t}"

Q = defaultdict(float)
penalty = 10.0  # constraint penalty weight

# ---- Constraint 1: Each operation assigned to exactly one time slot ----
for j, ops in jobs.items():
    for step, (m, duration) in enumerate(ops):
        # sum over t of x[j][m][t] == 1  =>  penalty (sum - 1)^2
        valid_times = range(max_time - duration + 1)
        for t in valid_times:
            v = var(j, m, t)
            Q[(v, v)] += -penalty  # linear: -1 * 2 * lambda (from expansion)
        for t1, t2 in itertools.combinations(valid_times, 2):
            v1, v2 = var(j, m, t1), var(j, m, t2)
            Q[(v1, v2)] += 2 * penalty  # quadratic cross terms

# ---- Constraint 2: Machine capacity (no overlap) ----
for m in range(n_machines):
    # Collect all jobs that use machine m
    machine_ops = [
        (j, step, t, duration)
        for j, ops in jobs.items()
        for step, (mach, duration) in enumerate(ops)
        if mach == m
        for t in range(max_time - duration + 1)
    ]
    for (j1, s1, t1, d1), (j2, s2, t2, d2) in itertools.combinations(machine_ops, 2):
        if j1 == j2:
            continue
        # Penalize if time intervals overlap
        if t1 < t2 + d2 and t2 < t1 + d1:
            v1, v2 = var(j1, m, t1), var(j2, m, t2)
            Q[(v1, v2)] += penalty

# ---- Objective: minimize makespan ----
# Add linear penalty for completion time of last operation per job
for j, ops in jobs.items():
    last_m, last_d = ops[-1]
    for t in range(max_time - last_d + 1):
        v = var(j, last_m, t)
        Q[(v, v)] += float(t + last_d) / max_time  # normalized completion time

# ---- Submit to D-Wave Advantage ----
bqm = dimod.BinaryQuadraticModel.from_qubo(Q)

sampler = EmbeddingComposite(DWaveSampler())
response = sampler.sample(
    bqm,
    num_reads=1000,        # number of annealing runs
    annealing_time=100,    # microseconds per anneal
    chain_strength=15.0,   # embedding chain coupling
    label="vw-jobshop-demo"
)

# Extract best valid solution
best_sample = None
best_makespan = float('inf')

for sample, energy in response.data(['sample', 'energy']):
    # Decode: find which time slot each operation was assigned
    schedule = defaultdict(dict)
    valid = True
    for j, ops in jobs.items():
        for step, (m, duration) in enumerate(ops):
            assigned = [t for t in range(max_time) if sample.get(var(j, m, t), 0) == 1]
            if len(assigned) != 1:
                valid = False
                break
            schedule[j][step] = (m, assigned[0], duration)
        if not valid:
            break
    if valid:
        makespan = max(
            schedule[j][len(ops)-1][1] + schedule[j][len(ops)-1][2]
            for j, ops in jobs.items()
        )
        if makespan < best_makespan:
            best_makespan = makespan
            best_sample = schedule

print(f"Best makespan found: {best_makespan}")

Embedding: Mapping the Problem onto Hardware

The D-Wave Advantage chip has a Pegasus graph topology: each qubit connects to at most 15 others. Real QUBO problems often require interactions between qubits that are not physically adjacent. The EmbeddingComposite class handles this automatically by chaining multiple physical qubits together to represent one logical variable.

For the Volkswagen factory scheduling problem with 50 stations and 30 jobs simultaneously on the line, the QUBO required roughly 4,500 logical variables and 22,000 quadratic terms. Embedding this onto Pegasus used about 3,800 physical qubits and 11,000 couplers, leaving room on the 5000+ qubit chip.

Results and Comparison

The Volkswagen team benchmarked D-Wave Advantage against:

• CPLEX branch-and-bound with 60-second time limit
• Simulated annealing (classical) with equivalent time budget
• Tabu search (production heuristic currently in use)

Results for 50-station instances averaged over 200 test cases:

Solver	Mean gap to optimum	Time to solution
CPLEX (60s)	1.1%	60 s
Simulated annealing	4.3%	5 s
Tabu search (production)	6.8%	2 s
D-Wave Advantage	2.7%	0.8 s

For instances with 80+ stations, CPLEX rarely found the optimum within 60 seconds, and D-Wave’s 2-3% gap was competitive with the best classical heuristics in the same time budget.

The most practically significant result: D-Wave could replan a disrupted schedule (one robot station goes offline) in under one second, while CPLEX-based replanning took 10-45 seconds. For a line running on 60-second cycle time, one second of schedule latency is acceptable; 45 seconds causes a cascade of delays.

Limitations and Honest Caveats

The results are promising but need careful interpretation:

• Solution quality depends heavily on penalty weight tuning. Different lambda values can dramatically change which solutions are found.
• Chain breaks (where physically chained qubits disagree) reduce solution quality. The team post-processed results by majority voting on chain values.
• D-Wave does not guarantee finding the global optimum. It is a heuristic, not an exact solver.
• Comparisons to classical methods depend on what time budget you allow the classical solver. Given unlimited time, classical exact methods will outperform annealing.
• The quantum speedup claim is disputed in academic literature. Several studies show well-tuned classical simulated annealing matches D-Wave on similar problems.

The practical case for D-Wave in this setting is not quantum speedup in theory, but a fast, parallelizable heuristic that integrates into real production systems.

Framework

D-Wave Ocean SDK is the primary development environment for quantum annealing:

pip install dwave-ocean-sdk
dwave setup  # configure API token and default solver

The Ocean SDK includes dimod (BQM construction), minorminer (embedding), and dwave-system (hardware access). Problems can also run on D-Wave’s classical hybrid solvers, which combine annealing with classical heuristics for very large instances.

Learn more: D-Wave Ocean SDK Reference