Quantum-classical hybrid optimisation for energy grids

Why quantum for energy grids

Power grid operation requires solving combinatorial optimisation problems daily: which generators to commit, at what output, over what time horizon, to meet demand at minimum cost while respecting ramp rates, minimum up/down times, and reserve requirements. This is the unit commitment (UC) problem — NP-hard in its general form.

Classical solvers (MILP via OR-Tools, Gurobi, CPLEX) handle moderate-size instances well, but the binary decision space grows exponentially: G generators over T time steps produce G×T binary variables. For large grids (hundreds of generators, 24–48 hour horizons), even commercial solvers hit practical limits or require aggressive relaxation.

Quantum computing offers a fundamentally different approach: encode the problem into a Hamiltonian whose ground state is the optimal commitment schedule, then use variational quantum algorithms (QAOA, VQE) to find low-energy states. QuantumGrid is a working implementation of this idea, benchmarked honestly against classical alternatives.

The unit commitment problem

Given a fleet of thermal generators, each with a capacity range [P_min, P_max], a marginal cost, startup/shutdown costs, and minimum up/down time constraints, the UC problem decides which generators to turn on (binary) at each hour to minimise total cost while meeting forecast demand.

The classical formulation is a Mixed-Integer Linear Program:

# Objective: minimise total cost
min  Σ_t Σ_g [ c_g · p_g,t + SU_g · y_g,t + SD_g · z_g,t ]

# Subject to:
Σ_g p_g,t ≥ D_t                          # demand balance
P_min_g · u_g,t ≤ p_g,t ≤ P_max_g · u_g,t  # capacity bounds
u_g,t - u_g,t-1 ≤ y_g,t                    # startup indicator
u_g,t-1 - u_g,t ≤ z_g,t                    # shutdown indicator
u_g,t ∈ {0, 1}                             # binary commitment

QuantumGrid solves this with Google OR-Tools (SCIP backend) as the classical reference. For our test instances (4 generators × 6 hours = 24 binary variables), MILP finds the proven optimum in under a second.

From MILP to QUBO to Ising

Quantum annealers and variational circuits don’t natively handle linear constraints and continuous variables. The standard pipeline converts MILP → QUBO → Ising:

• Step 1: Discretise or fix continuous power outputs (e.g., generators run at P_max when committed). This makes all variables binary.
• Step 2: Move constraints into the objective as penalty terms: H = H_cost + λ · H_penalty. Each violated constraint adds a quadratic penalty proportional to λ.
• Step 3: Express the resulting quadratic binary objective as a QUBO matrix Q where xTQx is to be minimised.
• Step 4: Map QUBO to Ising via the substitution x_i = (1 - z_i) / 2 where z_i ∈ {−1, +1}. This gives a Hamiltonian H = Σ_ij J_ij · z_i · z_j + Σ_i h_i · z_i + c.

def build_qubo(fleet, demand, horizon, penalty_lambda):
    n = len(fleet) * horizon  # number of binary variables
    Q = np.zeros((n, n))

    # Cost terms (diagonal)
    for g, gen in enumerate(fleet):
        for t in range(horizon):
            idx = g * horizon + t
            Q[idx, idx] += gen.marginal_cost * gen.p_max

    # Demand penalty (off-diagonal quadratic)
    for t in range(horizon):
        gen_indices = [g * horizon + t for g in range(len(fleet))]
        for i in gen_indices:
            Q[i, i] += penalty_lambda * (gen.p_max**2 - 2 * demand[t] * gen.p_max)
            for j in gen_indices:
                if i != j:
                    Q[i, j] += penalty_lambda * gen_i.p_max * gen_j.p_max

    return Q

The Ising mapping then converts Q into PennyLane observables: a sum of PauliZ ⊗ PauliZ coupling terms and single-qubit PauliZ fields.

QAOA: quantum approximate optimisation

The Quantum Approximate Optimisation Algorithm (Farhi et al., 2014) alternates between a cost unitary (encoding the problem Hamiltonian) and a mixer unitary (enabling exploration of the solution space). For depth p, the circuit has 2p variational parameters (γ₁, …, γ_p, β₁, …, β_p).

QuantumGrid implements QAOA using PennyLane’s ApproxTimeEvolution template for the cost unitary and RX gates for the X-mixer. The variational loop uses SciPy’s COBYLA optimiser (gradient-free, works well for noisy quantum objectives).

import pennylane as qml

dev = qml.device("default.qubit", wires=n_qubits)

@qml.qnode(dev)
def qaoa_circuit(gammas, betas, cost_hamiltonian):
    # Initial superposition
    for i in range(n_qubits):
        qml.Hadamard(wires=i)

    # p layers of cost + mixer
    for layer in range(p):
        qml.ApproxTimeEvolution(cost_hamiltonian, gammas[layer], 1)
        for i in range(n_qubits):
            qml.RX(2 * betas[layer], wires=i)

    return qml.expval(cost_hamiltonian)

# Optimise with COBYLA
from scipy.optimize import minimize
result = minimize(
    lambda params: qaoa_circuit(params[:p], params[p:], H_cost),
    x0=np.random.uniform(0, np.pi, 2 * p),
    method="COBYLA",
    options={"maxiter": 200}
)

With p=1 on a 24-qubit instance (4 generators × 6 hours), QAOA completes in ~100 seconds on CPU using the statevector simulator. Solution quality varies with random initialisation — multiple restarts are essential.

VQE: variational quantum eigensolver

VQE (Peruzzo et al., 2014) uses a hardware-efficient ansatz instead of the problem-specific QAOA structure. QuantumGrid implements an R_y/R_z + CNOT ladder ansatz with configurable depth. Each qubit gets R_y and R_z rotations, followed by a linear chain of CNOTs for entanglement.

@qml.qnode(dev)
def vqe_circuit(params, cost_hamiltonian):
    n_layers = params.shape[0] // (2 * n_qubits)
    idx = 0
    for layer in range(n_layers):
        for i in range(n_qubits):
            qml.RY(params[idx], wires=i); idx += 1
            qml.RZ(params[idx], wires=i); idx += 1
        for i in range(n_qubits - 1):
            qml.CNOT(wires=[i, i + 1])
    return qml.expval(cost_hamiltonian)

VQE is more expressive than QAOA but scales poorly: for 24 qubits with 1 layer, there are 48 parameters (vs 2 for QAOA p=1). The optimisation landscape is harder, and barren plateaus become a concern for larger instances.

Classical baselines

Honest benchmarking requires strong classical baselines. QuantumGrid implements three:

For the small instances where quantum solvers can run (6–24 qubits), MILP always finds the optimal solution. The point of quantum isn’t to beat MILP on 24 variables — it’s to demonstrate the formulation pipeline and study how variational algorithms behave as the problem scales toward regimes where classical solvers struggle.

Penalty tuning

The QUBO formulation moves constraints into the objective via penalty weight λ. If λ is too small, the solver finds low-cost but infeasible solutions (generators don’t meet demand). If λ is too large, the penalty dominates and the solver can’t distinguish between feasible solutions of different cost — it effectively ignores the cost function.

QuantumGrid implements binary search for λ: start with a range [λ_min, λ_max], solve at the midpoint, check feasibility, and narrow the range. The sensitivity sweep then evaluates solution quality across a grid of λ values to find the stable region.

# Penalty tuning via binary search
lambda_min, lambda_max = 0.1, 100.0
for _ in range(20):
    lam = (lambda_min + lambda_max) / 2
    Q = build_qubo(fleet, demand, horizon, penalty_lambda=lam)
    solution = solve_qaoa(Q, p=1)
    if is_feasible(solution, demand):
        lambda_max = lam       # try smaller penalty
    else:
        lambda_min = lam       # need stronger enforcement
best_lambda = (lambda_min + lambda_max) / 2

Scaling analysis

How do quantum and classical solvers scale as the problem grows? QuantumGrid runs a systematic sweep across problem sizes:

The statevector simulator’s memory grows as 2ⁿ: 24 qubits requires 128 MB, 30 qubits would need 8 GB. Beyond ~26 qubits, you need either a tensor-network simulator or real quantum hardware. MILP, by contrast, routinely handles hundreds of binary variables.

The honest conclusion: on instances that fit in a statevector simulator, classical solvers dominate. The potential quantum advantage lies at scales where the exponential blowup of branch-and-bound meets the (hoped-for) polynomial scaling of quantum optimisation — but we’re not there yet.

Architecture diagram

The full QuantumGrid pipeline from raw data to solved commitment schedule:

Key takeaways

• The QUBO pipeline is the real contribution. Converting a real-world MILP into a quantum-compatible formulation is non-trivial. The encoding, penalty calibration, and Ising mapping are where most of the engineering lives — the quantum circuit is the easy part.
• Penalty tuning is critical. A factor-of-2 error in λ takes you from feasible solutions to garbage. Binary search with feasibility checks is the minimum viable approach; Lagrangian relaxation would be more principled.
• QAOA scales better than VQE for structured problems. QAOA with p=1 has only 2 parameters regardless of qubit count. VQE’s hardware-efficient ansatz grows linearly with qubits, making optimisation harder and more prone to barren plateaus.
• Classical solvers are still faster. For any problem size that fits on a statevector simulator (~26 qubits), MILP solves it optimally in milliseconds. The quantum advantage regime, if it exists, is at scales we can’t yet simulate classically.
• Honest benchmarking matters. It’s easy to publish a QAOA paper that shows “competitive results” by choosing problem instances where all solvers do well. QuantumGrid benchmarks all 5 solvers on the same instances with the same metrics.
• Energy grids are a natural quantum target. The unit commitment problem has the right structure: binary decisions, quadratic objectives, constraints that map to penalty terms. Whether quantum hardware will deliver a practical speedup depends on error rates and qubit counts that don’t yet exist — but the formulation pipeline is ready.

Full source, 130+ tests across 10 releases, and benchmark pipeline at github.com/ajliouat/quantumgrid.