Siemens AG: Quantum Optimization for Industrial Job Shop Scheduling

The Problem

Siemens manufactures complex industrial systems: trains, gas turbines, factory automation equipment. A production schedule involves thousands of tasks, hundreds of machines, and dependency constraints spanning weeks. Finding the optimal schedule minimizes makespan (total production time). Job shop scheduling is NP-hard, and industrial instances far exceed the scale of classical exact solvers.

Siemens’ quantum team asked whether quantum solvers could contribute on the bottleneck subproblems where classical search is most expensive.

QUBO Formulation

Binary variable x[j, op, t] = 1 means operation op of job j starts at time t. Three penalty terms encode the constraints: each operation starts exactly once, operations within a job respect precedence order, and no two operations share a machine at the same time.

import dimod
import numpy as np

# 3-job, 3-machine instance: each job is [(machine_id, duration), ...]
jobs = [
    [(0, 3), (1, 2), (2, 2)],
    [(1, 2), (0, 3), (2, 1)],
    [(2, 2), (1, 1), (0, 2)],
]
n_machines, n_time_slots = 3, 10

def build_job_shop_qubo(jobs, n_time_slots, A=10.0, B=8.0, C=6.0):
    """QUBO penalties: A=one-start, B=precedence, C=no-overlap."""
    Q = {}
    def var(j, op, t):
        return (sum(len(jobs[jj]) for jj in range(j)) + op) * n_time_slots + t

    # Constraint 1: each operation starts exactly once
    for j, job in enumerate(jobs):
        for op in range(len(job)):
            for t in range(n_time_slots):
                v = var(j, op, t)
                Q[(v, v)] = Q.get((v, v), 0) + A * (1 - 2)
                for t2 in range(t + 1, n_time_slots):
                    v2 = var(j, op, t2)
                    Q[(min(v, v2), max(v, v2))] = Q.get((min(v, v2), max(v, v2)), 0) + 2 * A

    # Constraint 2: precedence within each job
    for j, job in enumerate(jobs):
        for op in range(len(job) - 1):
            for ts in range(n_time_slots):
                for tn in range(min(ts + job[op][1], n_time_slots)):
                    v1, v2 = var(j, op, ts), var(j, op + 1, tn)
                    Q[(min(v1, v2), max(v1, v2))] = Q.get((min(v1, v2), max(v1, v2)), 0) + B

    for j1, job1 in enumerate(jobs):
        for op1, (m1, d1) in enumerate(job1):
            for j2, job2 in enumerate(jobs):
                if j2 <= j1:
                    continue
                for op2, (m2, d2) in enumerate(job2):
                    if m1 != m2:
                        continue
                    for t1 in range(n_time_slots):
                        for t2 in range(n_time_slots):
                            if t1 < t2 + d2 and t2 < t1 + d1:
                                v1, v2 = var(j1, op1, t1), var(j2, op2, t2)
                                Q[(min(v1, v2), max(v1, v2))] = (
                                    Q.get((min(v1, v2), max(v1, v2)), 0) + C)

    return dimod.BinaryQuadraticModel.from_qubo(Q), var

bqm, var_fn = build_job_shop_qubo(jobs, n_time_slots)
print(f"QUBO: {len(bqm.variables)} variables, {len(bqm.quadratic)} interactions")

Solving with D-Wave Hybrid and Decoding

from dwave.system import LeapHybridSampler

result = LeapHybridSampler().sample(bqm, time_limit=10, label="siemens-job-shop")
sample = result.first.sample

# Decode start times
schedule = {j: {op: t for op in range(len(job)) for t in range(n_time_slots)
                if sample.get(var_fn(j, op, t), 0) == 1}
            for j, job in enumerate(jobs)}
makespan = max(schedule[j].get(op, 0) + jobs[j][op][1]
               for j, job in enumerate(jobs) for op in range(len(job)))
print(f"Energy: {result.first.energy:.2f}  Makespan: {makespan}")
for j, job in enumerate(jobs):
    for op, (machine, dur) in enumerate(job):
        print(f"  Job {j} op {op}: machine {machine} start={schedule[j].get(op, '?')}")

Classical Decomposition Architecture

For large instances beyond current qubit capacity, Siemens used classical decomposition: find the bottleneck machine (highest total load), extract the scheduling subproblem for just those operations, send it to the quantum solver, and assemble the full schedule classically around the quantum result.

machine_loads = {m: 0 for m in range(n_machines)}
for job in jobs:
    for machine, duration in job:
        machine_loads[machine] += duration

bottleneck = max(machine_loads, key=machine_loads.get)
bottleneck_ops = [(j, op, dur) for j, job in enumerate(jobs)
                  for op, (machine, dur) in enumerate(job) if machine == bottleneck]
print(f"Bottleneck machine {bottleneck} (load={machine_loads[bottleneck]})")
print(f"Subproblem: {len(bottleneck_ops)} ops -> quantum solver")
print("Remaining operations -> classical scheduling around quantum result")

Results

Siemens compared D-Wave hybrid, QAOA (p=1,2,3), and Gurobi on 10-50 job instances:

• 10 jobs: QAOA at p=2 within 5% of Gurobi optimal; D-Wave hybrid matched Gurobi exactly
• 20-50 jobs: QAOA infeasible (too many qubits); D-Wave hybrid competitive in quality and faster on dense instances
• Hybrid decomposition scaled to 200-job factory instances with quantum solvers handling 10-15 job subproblems

Full end-to-end advantage on factory scheduling requires 10,000+ logical qubits. Siemens estimates that milestone is a decade away, but the hybrid decomposition methodology is applicable on current hardware today.

Learn more: D-Wave Ocean SDK Reference