ExxonMobil: Quantum Optimisation for Shipping Route Planning

The Problem

Maritime shipping is one of the most fuel-intensive industries on Earth. Fuel accounts for roughly 50% of total operating costs for a large LNG tanker, and ExxonMobil operates a global fleet spanning dozens of vessels, hundreds of port calls, and tight contractual delivery windows. Optimising routes across this network is not a simple shortest-path problem.

The underlying mathematical structure is a variant of the Vehicle Routing Problem (VRP): assign a fleet of vessels to delivery routes such that total fuel consumption is minimised, subject to constraints on cargo capacity, port berthing windows, vessel availability, and weather routing corridors. VRP is NP-hard. For a fleet of 20 tankers serving 40 ports with varying time windows, the solution space contains more candidate assignments than can be exhaustively searched by any classical computer in a reasonable timeframe.

Classical solvers (CPLEX, Gurobi, simulated annealing) attack VRP through branch-and-bound, branch-and-cut, or metaheuristic search. They work well up to moderate problem sizes, but as fleet size and port network density increase, runtime grows quickly and solution quality degrades. ExxonMobil’s question to IBM Quantum was direct: can quantum algorithms find better solutions, or find equivalent solutions faster, as qubit counts and gate fidelities improve?

QUBO Formulation

Before a problem reaches a quantum computer, it must be cast as a Quadratic Unconstrained Binary Optimisation (QUBO) problem. Each binary variable represents a decision: x_{i,r} = 1 means vessel i is assigned to route r.

The objective encodes fuel cost as a sum of linear and quadratic terms. Constraints (each port visited exactly once, vessel capacity limits, time window feasibility) enter as penalty terms with a Lagrange multiplier lambda large enough to make infeasible solutions energetically unfavourable.

A 10-vessel, 20-port instance with discretised time windows maps to roughly 200-400 binary variables after encoding. The QUBO matrix Q captures all pairwise interactions between these variables.

From QUBO, the step to an Ising Hamiltonian is direct. Substituting x_i = (1 - Z_i) / 2, where Z_i is a Pauli-Z operator acting on qubit i, converts the binary optimisation into a sum of Z and ZZ Pauli terms acting on a register of qubits. This cost Hamiltonian H_C is then the target for QAOA.

QAOA Circuit Construction in Qiskit

QAOA alternates between a cost unitary exp(-igammaH_C) and a mixer unitary exp(-ibetaH_X), where H_X is a sum of Pauli-X operators. The number of alternating layers is the circuit depth p. Larger p gives a better approximation but requires more gates and more coherence time.

from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel
from qiskit_ibm_runtime import QiskitRuntimeService
import numpy as np
from scipy.optimize import minimize

# Problem: 10 nodes (vessels + ports), small benchmark instance
n_qubits = 10
p = 2  # QAOA layers; p=1 tested first, p=2 showed improvement

# QUBO coefficients (diagonal = linear, off-diagonal = quadratic)
# In practice these come from the routing cost matrix and penalty terms
Q = np.random.randn(n_qubits, n_qubits)
Q = (Q + Q.T) / 2  # symmetrise

def build_qaoa_circuit(gamma, beta, Q, n, p):
    qc = QuantumCircuit(n, n)

    # Initialise in equal superposition
    qc.h(range(n))

    for layer in range(p):
        # Cost unitary: ZZ terms from off-diagonal Q, Z terms from diagonal
        for i in range(n):
            for j in range(i + 1, n):
                if abs(Q[i, j]) > 1e-8:
                    qc.rzz(2 * gamma[layer] * Q[i, j], i, j)
            if abs(Q[i, i]) > 1e-8:
                qc.rz(2 * gamma[layer] * Q[i, i], i)

        # Mixer unitary: transverse field X on every qubit
        for i in range(n):
            qc.rx(2 * beta[layer], i)

    qc.measure(range(n), range(n))
    return qc

def evaluate_qubo(bitstring, Q):
    x = np.array([int(b) for b in bitstring])
    return float(x @ Q @ x)

def qaoa_objective(params, Q, n, p, shots=2048):
    gamma = params[:p]
    beta = params[p:]
    qc = build_qaoa_circuit(gamma, beta, Q, n, p)
    sim = AerSimulator()
    counts = sim.run(qc, shots=shots).result().get_counts()
    # Expected energy over sampled bitstrings
    total_energy = sum(
        evaluate_qubo(bs, Q) * count
        for bs, count in counts.items()
    )
    return total_energy / shots

# Optimise variational parameters with COBYLA
x0 = np.random.uniform(0, np.pi, 2 * p)
result = minimize(qaoa_objective, x0, args=(Q, n_qubits, p),
                  method='COBYLA', options={'maxiter': 500})
print(f"Optimised QAOA energy: {result.fun:.4f}")

The p=1 circuit depth requires approximately 45 two-qubit gates for a 10-qubit fully-connected QUBO instance. The p=2 circuit doubles this. On IBM Quantum Eagle hardware, the two-qubit gate error rate sits around 0.1-0.3%, meaning a 90-gate circuit accumulates meaningful noise before measurement.

Benchmarking Against Classical Solvers

ExxonMobil and IBM benchmarked QAOA against two classical baselines on instances of 10-20 nodes:

• CPLEX: IBM’s commercial integer programming solver, which finds provably optimal solutions for small instances
• Simulated Annealing: a metaheuristic that provides high-quality heuristic solutions quickly

On 10-node instances, QAOA with p=2 found solutions within 2-5% of the CPLEX optimum on the noiseless simulator. On real Eagle hardware, noise degraded solution quality such that QAOA solutions were 8-15% above optimum, comparable to or slightly worse than simulated annealing run for the same wall-clock time.

On 20-node instances (still tiny relative to ExxonMobil’s real fleet), the gap widened. Hardware noise accumulates with circuit depth, and QAOA at p=2 on 20 qubits involves enough gates that error rates become the dominant factor in solution quality.

The honest summary: QAOA is not superior to classical methods on any practically sized routing problem using 2023 hardware. The results confirm theoretical expectations but do not demonstrate practical advantage.

Error Mitigation Techniques

The team applied several noise mitigation strategies to improve hardware results:

Zero-noise extrapolation (ZNE): Run the circuit at artificially increased noise levels (by repeating gates) and extrapolate to the zero-noise limit. This adds circuit executions but reduces systematic bias in expectation values.

Measurement error mitigation: Characterise readout errors on each qubit with a calibration matrix and invert it to correct measured bitstring distributions.

Dynamical decoupling: Insert sequences of X gates during idle periods to suppress decoherence from environmental noise.

Even with all three applied, hardware results on 20-qubit instances showed 10-12% above optimum, compared to 3-4% on the noiseless simulator. The gap closes noticeably with error mitigation but does not vanish.

Resource Estimation for Fault-Tolerant Advantage

IBM and ExxonMobil’s joint analysis estimated what fault-tolerant hardware would be required to achieve meaningful quantum advantage on a real fleet routing problem.

A realistic ExxonMobil routing instance might involve 50 vessels and 200 ports with time window constraints, encoding to roughly 10,000 binary variables after problem reduction. Running QAOA at p=10-15 (the depth range where theory suggests advantage over simulated annealing for this problem class) on 10,000 logical qubits requires:

• Physical qubit count: 1-10 million (depending on error correction code and target logical error rate)
• Gate counts: tens of millions of logical operations
• Runtime per circuit evaluation: minutes to hours on projected fault-tolerant hardware

Current estimates from IBM’s roadmap place fault-tolerant devices with thousands of logical qubits in the 2030-2033 timeframe. Practical routing advantage, requiring tens of thousands of logical qubits and low logical error rates, is a decade or more away.

IBM and ExxonMobil’s 10-Year Roadmap

The collaboration laid out a phased plan:

Phase 1 (2021-2024): Algorithm development and NISQ benchmarking. Identify which problem substructures are most amenable to near-term quantum approaches. Develop error mitigation workflows. This phase is complete.

Phase 2 (2025-2028): Early fault-tolerant experiments. As IBM releases devices with tens to hundreds of logical qubits (Condor, Heron, and successor chips), run progressively larger routing instances. Develop quantum-classical hybrid pipelines that use quantum subroutines within a classical VRP solver framework.

Phase 3 (2029+): Full fault-tolerant advantage. If hardware milestones are met, quantum algorithms may deliver demonstrable advantage on real fleet routing. The target metric is fuel savings: even a 1% improvement across ExxonMobil’s fleet represents hundreds of millions of dollars annually.

What This Means for Practitioners

The ExxonMobil work is a serious, methodologically rigorous demonstration that QAOA can be applied to industrial logistics problems. It is also an honest acknowledgement that current hardware is not ready.

For quantum computing students, this case study illustrates several important lessons. First, the QUBO and Ising Hamiltonian formulation is the standard entry point for combinatorial optimisation on quantum hardware. Second, circuit depth is the primary enemy of solution quality on NISQ devices, and problem-specific compilation matters. Third, benchmarking against mature classical solvers rather than against naive baselines is what separates credible research from hype.

The routing problem will eventually be a strong candidate for quantum advantage, because the problem scales hard and because even marginal improvements in solution quality translate to large real-world savings. But “eventually” is the operative word.

Learn more: Qiskit Reference