SAP Quantum Supply Chain Optimization for Manufacturing

The Problem

SAP’s supply chain planning software, SAP Integrated Business Planning (IBP), serves manufacturers managing hundreds of distribution facilities and thousands of SKUs. At that scale, multi-echelon inventory optimization (deciding how much safety stock to hold at each node in a supply network to meet service level targets at minimum cost) is an NP-hard combinatorial problem. So is the related supplier selection problem: choosing which suppliers to certify for each component given cost, lead time, and quality constraints across a global bill of materials.

Classical approaches use mixed-integer programming (MIP) solvers such as CPLEX or Gurobi, which work well for problems up to a few hundred variables but slow dramatically as problem size grows. For the largest SAP customers, the planning horizon combined with the number of nodes produces MIP instances that take hours to solve, constraining how frequently plans can be refreshed. SAP began evaluating quantum and quantum-inspired solvers as a path to better solutions at acceptable runtimes.

QUBO Formulation

The core of both inventory optimization and supplier selection is reformulating the problem as a quadratic unconstrained binary optimization (QUBO), the native input format for both quantum annealers and QAOA circuits. For the supplier selection problem, each binary variable x_ij encodes whether supplier i is certified for component j. Cost terms penalize total certification cost; constraint terms penalize configurations that leave any component without a qualified supplier.

Integer variables (order quantities, inventory levels) require binary encoding: a variable with range [0, 63] requires 6 binary variables. For multi-echelon inventory, this encoding multiplies the effective problem size, making the gap between QUBO variable count and the original MIP variable count roughly 4-8x.

from qiskit_optimization import QuadraticProgram
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
import numpy as np

# Simplified 12-node supply network: 4 suppliers, 4 distribution centers, 4 customers
n_suppliers = 4
n_dcs = 4
n_customers = 4

# Cost matrix: cost[i][j] = unit transport cost from node i to node j
np.random.seed(42)
transport_cost = np.random.uniform(1.0, 10.0, (n_suppliers, n_dcs))
holding_cost = np.random.uniform(0.5, 3.0, n_dcs)
demand = np.random.uniform(100, 500, n_customers)

# Build quadratic program for supplier selection subproblem
qp = QuadraticProgram("SupplierSelection")

# Binary variable x[i][j]: supplier i assigned to serve DC j
for i in range(n_suppliers):
    for j in range(n_dcs):
        qp.binary_var(f"x_{i}_{j}")

# Objective: minimize total transport cost
linear = {}
for i in range(n_suppliers):
    for j in range(n_dcs):
        linear[f"x_{i}_{j}"] = transport_cost[i, j]
qp.minimize(linear=linear)

# Constraint: each DC must be served by exactly one supplier
for j in range(n_dcs):
    constraint = {f"x_{i}_{j}": 1 for i in range(n_suppliers)}
    qp.linear_constraint(linear=constraint, sense="=", rhs=1, name=f"dc_{j}_covered")

# Convert to QUBO (penalty method for equality constraints)
converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)

print(f"Original variables: {qp.get_num_vars()}")
print(f"QUBO variables: {qubo.get_num_vars()}")
print(f"QUBO linear terms: {len(qubo.objective.coefficients.to_dict())}")

# Run QAOA with p=2 layers
sampler = Sampler()
qaoa = QAOA(sampler=sampler, optimizer=COBYLA(maxiter=200), reps=2)
optimizer = MinimumEigenOptimizer(qaoa)

result = optimizer.solve(qubo)
print(f"QAOA objective value: {result.fvalue:.2f}")
print(f"Selected assignment: {result.x}")

D-Wave Hybrid Solver for Larger Instances

For the 50-node test case (20 suppliers, 15 distribution centers, 15 customers), QAOA on IBM Falcon was too shallow to produce competitive solutions; the QUBO had 300 binary variables after integer encoding, far beyond what current NISQ QAOA handles well. SAP switched to D-Wave’s hybrid solver, which runs a classical branch-and-bound search guided by D-Wave Advantage’s quantum annealer for subproblem evaluation.

The D-Wave Advantage processor has over 5000 qubits, but with limited connectivity (the Pegasus graph topology). Embedding the 300-variable dense QUBO onto the sparse Pegasus graph required chain lengths of 3-5 physical qubits per logical variable, effectively reducing the usable logical qubit count. D-Wave’s hybrid BQM (binary quadratic model) solver handled this embedding automatically.

import dimod
import dwave.system as dwave_sys
from dwave.system import LeapHybridSampler

# Construct BQM from the 50-node supply chain QUBO
# (linear and quadratic coefficients from the full problem formulation)
bqm = dimod.BinaryQuadraticModel(vartype="BINARY")

# Add variables and interactions (simplified illustration)
n_vars = 300
for i in range(n_vars):
    bqm.add_variable(i, np.random.uniform(-1, 1))  # linear bias

for i in range(n_vars):
    for j in range(i + 1, min(i + 5, n_vars)):     # sparse interactions
        bqm.add_interaction(i, j, np.random.uniform(-0.5, 0.5))

# D-Wave Leap hybrid solver: classical + quantum annealing
sampler = LeapHybridSampler()
response = sampler.sample(bqm, time_limit=60)  # 60-second solve budget

best = response.first
print(f"Best energy found: {best.energy:.4f}")
print(f"Number of variables: {len(best.sample)}")

SAP IBP Integration Roadmap

SAP’s architecture team evaluated three integration patterns for quantum solvers within SAP IBP: direct API calls to cloud quantum backends during the planning run (latency-sensitive, requires fast queue access), pre-computation of quantum-optimized parameters used to warm-start classical MIP solvers (the current preferred approach), and quantum-in-the-loop optimization for real-time re-planning triggered by supply disruptions.

The 15% cost reduction observed with the D-Wave hybrid solver on the 50-node instance was reproducible across five independent runs, providing confidence that the result reflects genuine solution quality improvement rather than noise. On the 12-node instance where QAOA results could be compared against exact solutions, QAOA matched the optimum in 9 of 10 runs with p=2 layers and COBYLA optimization.

SAP published findings through the IBM Quantum Network research program and filed a patent application on the QUBO formulation method for multi-echelon inventory with service level constraints.

Learn more: D-Wave Ocean Reference | Qiskit Optimization Reference