Quantum computing algorithms for two-stage stochastic optimization, with a focus on the Unit Commitment (UC) problem in power systems. The algorithms combine Discrete Quantum Annealing (DQA) with Quantum Amplitude Estimation (QAE) to compute expected-value objective functions over a probability distribution of wind-power scenarios, which is the key subroutine needed for first-stage hedged decision-making.
Based on: arXiv 2402.15029 — "Quantum algorithms for the two-stage stochastic unit commitment problem"
- Problem description
- Repository structure
- Module reference —
qiskit_impl/ - Module reference —
QuantumExpectedValueFunctionProject/ - Installation and environment
- Quick-start usage
- Notebooks
- GPU porting strategy (branch
gpuOptim) - Fault-tolerant resource estimation
- Development conventions
- References
The Unit Commitment problem decides which gas generators to commit before uncertain wind power is realised (first-stage decision), then dispatches all generators optimally for each realised scenario (second-stage recourse).
Mathematically:
-
First-stage: binary gas commitment
$x \in {0,1}^{n_g}$ -
Second-stage recourse: real-valued
$r(\xi) \ge 0$ to cover unmet demand at cost$c_r$ -
Uncertainty: wind scenario
$\xi \sim p(\xi)$ (encoded as a quantum PDF register)
The quantum algorithm evaluates
quantum_stochastic_programming/
│
├── qiskit_impl/ # Modular Qiskit 2.x implementation (primary)
│ ├── binary_optimizer.py # Core DQA+QAE solver (BinaryNestedOptimizer)
│ ├── qae.py # Standalone QAE circuit builder (QAE_Optimizer)
│ ├── ExpValFun_functions.py # Pure-function circuit primitives
│ ├── dist_prep.py # PDF encoding circuits
│ ├── resource_estimator.py # Fault-tolerant resource estimation
│ └── qae_example.ipynb # Example notebook
│
├── QuantumExpectedValueFunctionProject/ # Original implementation (C. Rotello)
│ ├── optimizer_utils.py # VariableRegister + PowerSystem_1Bus
│ ├── dense_optimizer.py # Optimizer_Dense — scenario-superposition solver
│ ├── expanded_optimizer.py # Optimizer_Expanded — scenario-expanded solver
│ ├── binary_optimizer.py # Legacy BinaryNestedOptimizer
│ ├── experiment_convergence.py # Experiment runner script
│ ├── test_dense_optimizer.py # Unit tests for dense_optimizer
│ ├── test_expanded_optimizer.py # Unit tests for expanded_optimizer
│ ├── test_optimizer_utils.py # Unit tests for optimizer_utils
│ ├── four_binary_turbines_uniform_distribution_tutorial.ipynb # Tutorial
│ ├── ExpValFun_QAOA_QAE_computations.ipynb # Computation notebook
│ ├── ExpValFun_QAOA_QAE_figures.ipynb # Figure generation
│ ├── simple_UC_model.ipynb # Minimal UC example
│ └── SED_classical_solution.ipynb # Classical reference solution
│
├── cuda-q_mac_install_notes.md # CUDA-Q environment notes
├── twoStageStocOptQ.code-workspace # VS Code multi-root workspace
└── README.md
Two implementations exist because the code evolved:
qiskit_impl/ |
QuantumExpectedValueFunctionProject/ |
|
|---|---|---|
| API style | Functional primitives + args dicts |
Class-based OOP |
| Encoding | Binary only | Binary or Unary |
| QAE integration | QFTGate-based canonical QAE | Multiple annealing solve strategies |
| Qiskit version | 2.x (qiskit_aer separated) | Originally 0.x/1.x (needs import fixes) |
| Status | Active development | Legacy / reference |
The main entry point for the DQA+QAE pipeline.
Wraps a DQA state-preparation circuit execute_qae().
| Argument | Type | Description |
|---|---|---|
uopt |
QuantumCircuit |
State-preparation unitary (DQA output) |
oracle |
QuantumCircuit |
Phase-marking oracle |
n |
int |
Number of wind qubits |
The core solver. Encodes the power system as a quantum circuit and drives the full DQA → QAE pipeline.
bno = BinaryNestedOptimizer(
gas_costs = [0.4], # first-stage generator costs
wind_costs = [0.05, 0.10], # wind turbine recourse cost coefficients
recourse_cost = 1.0, # penalty for unmet demand
pdf = {(0,0):0.25, (0,1):0.25, (1,0):0.25, (1,1):0.25},
demand = 1, # total wind demand (sum of y_j)
is_uniform = False, # True → use uniform distribution shortcut
)Circuit-building methods:
| Method | Returns | Description |
|---|---|---|
dicke_state_initialize(weight) |
QuantumCircuit |
Prepares weight excitations (feasible initial state) |
pdf_initialize() |
QuantumCircuit |
Encodes |
cost_operator(amplitude, constraint_amplitude, norm) |
QuantumCircuit |
Phase operator |
cost_scenario_operator(scenario, amplitude, constraint_amplitude, norm) |
QuantumCircuit |
Single-scenario phase operator |
demand_constraint_preserving_mixer(amplitude) |
QuantumCircuit |
XY-mixer |
adiabatic_evolution_circuit(wind_demand, time, time_steps, norm) |
QuantumCircuit |
Full DQA circuit |
adiabatic_evolution_scenario_circuit(wind_demand, time, time_steps, scenario, norm) |
QuantumCircuit |
Single-scenario DQA circuit |
implemented_qae(op, oracle, op_inv, oracle_inv, m, norm) |
QuantumCircuit |
Canonical QAE circuit (Brassard 2002, Fig. 2 of paper) |
canonical_qae(Uopt, oracle, m, c, norm) |
QuantumCircuit |
Convenience wrapper for canonical QAE |
Oracle methods (mark the good subspace for QAE):
| Method | Description |
|---|---|
exact_oracle(ydemand, norm, inverse) |
Exact phase-kickback oracle |
single_oracle_sin_inconstraint(c, norm, inverse) |
|
single_oracle_asin_inconstraint(norm, inverse) |
|
grover_operator_sin_inconstraint(power, c, norm) |
Grover operator |
grover_reflections_sin_inconstraint(m, c, norm) |
|
Execution methods:
| Method | Description |
|---|---|
execute_optimizer(qc, num_meas) |
Run DQA circuit on statevector_simulator, return measurement counts |
execute_optimizer_oracle(qc, num_meas) |
Run oracle circuit on aer_simulator_statevector
|
execute_qae(qc, m, num_meas) |
Run canonical QAE circuit with |
Classical utilities:
| Method | Description |
|---|---|
gen_all_wind_outputs() |
Enumerate all feasible wind assignments |
wind_scenario_cost(wind_output, scenario, wind_demand) |
Evaluate cost for one assignment + scenario |
brute_force_wind_demand_expectation_values() |
Classical brute-force expected value for all demands (reference) |
brute_force_energy_surface() |
Full energy landscape (for small instances) |
prep_gs(wind_demand) |
Prepare the exact ground-state statevector (for testing) |
process_expectation_value_optimizer(wind_demand, counts) |
Convert measurement counts to expected cost |
Standalone QAE circuit builder, separated from the optimizer logic.
Assembles the canonical QAE circuit from an args dict.
args = {
'n_y' : 4, # number of wind qubits
'm' : 3, # number of evaluation (ancilla) qubits
'pdf' : pdf_dict, # scenario probability distribution
'y_reg' : [0,1,2,3], # qubit indices for wind register
'pdf_reg': [4,5,6,7], # qubit indices for PDF register
'Theta' : theta, # QAE target angle
'w_d' : wind_demand, # Hamming weight constraint
'norm' : 10.0, # cost normalization
'uniform': False, # use uniform PDF shortcut
'gateset': 'default',
'cost_operator_circuit' : cost_op_fn, # callable -> QuantumCircuit
'mixer_operator_circuit' : mixer_op_fn, # callable -> QuantumCircuit
}
qae_opt = QAE_Optimizer(args)
qc = qae_opt.compile_qae_circuit() # full QAE circuit
qc = qae_opt.implemented_qae(op, oracle, op_inv, oracle_inv) # lower-levelPure functions returning QuantumCircuit objects. No class state — designed for use as callables passed via args dicts.
| Function | Signature | Description |
|---|---|---|
demand_constraint_preserving_mixer |
(amplitude, args) |
XY mixer; preserves w_d
|
cost_operator |
(amplitude, args) |
Full stochastic cost phase operator |
cost_scenario_operator |
(amplitude, args) |
Single-scenario cost phase operator |
single_oracle_sin_inconstraint |
(args, inverse=False) |
|
dicke_state_circuit |
(args) |
Dicke state preparation |
alternating_operator_ansatz |
(args) |
Full DQA alternating-operator ansatz for one Trotter step |
Circuits that encode probability distributions into quantum amplitudes.
| Function | Description |
|---|---|
pdf_initialize(args) |
Encode arbitrary discrete PDF |
make_pdf_uniform(n) |
Uniform distribution over |
make_pdf_skewnormal(n) |
Skew-normal distribution, |
make_normal_distribution_circuit(args) |
Gaussian amplitude encoding |
make_variational_distribution_circuit(args) |
Variational (trainable) amplitude encoding |
Fault-tolerant resource estimation using pytket (gate counting) and qualtran (Beverland et al. surface-code model).
args = {
'num_algo_qubits' : 10,
'circuit' : pytket_circuit,
'num_meas' : 1000,
'error_budget' : 1e-3,
'architecture' : 'superconducting', # or 'ion', 'majorana'
'optimistic' : False,
# Magic state factory params:
'R' : 2, # distillation rounds
'Q_0' : 1e-3, # physical qubit error rate
'P_r' : 1e-6, # target logical Clifford error
'F_copies': 1,
}
re = Resource_Estimator(args)
result = re.bev_resource_estimator()
# result keys: num_physical_qubits, num_logical_qubits, T_count, T_depth,
# clifford_count, rotation_count, wall_clock_time, ...Module-level helpers:
| Function | Description |
|---|---|
get_circuit_metrics(circuit) |
Extract gate counts (T, Clifford, rotation) from a pytket Circuit |
get_factory_params(R, Q_0, P_r, d, c, M, tau_d, C, error_budget) |
Compute magic-state factory parameters |
qasm_to_clifford_and_t(qc, ...) |
Convert Qiskit QuantumCircuit → Clifford+T decomposition via Solovay-Kitaev |
Foundation utilities shared by all optimizers.
Stores an integer variable as a qubit register. Handles binary and unary encodings, and generates Hamiltonian sub-circuits for arithmetic constraints.
reg = VariableRegister(max_value=3, encoding='binary') # 2-qubit registerKey methods:
| Method | Description |
|---|---|
numberOperator(amplitude) |
Phase operator |
squaredOperator(amplitude) |
Phase operator |
productOperator(other, amplitude) |
Two-register product |
swapOperator(other, amplitude) |
SWAP-network mixing between two registers |
lessThanOperator(other, amplitude) |
Phase if self < other
|
lessThanValue(value, amplitude) |
Phase if self < value
|
getValue(bstr) |
Decode measurement bitstring to integer |
setValue(integer) |
Encode integer as circuit preparation |
Stores the power-system problem specification.
system = PowerSystem_1Bus(
gas_costs = [0.4],
wind_costs = [0.05, 0.10],
decision_levels = 4,
undersatisfied_cost = 1.0,
demand = 2,
pdf = {(0,1): 0.5, (1,0): 0.5},
normalization = None,
)Key methods:
| Method | Description |
|---|---|
cobylaSolve(discrete, gas_values) |
Classical COBYLA solver (reference) |
getFirstStageCosts(gas_decisions) |
Evaluate first-stage cost for a gas commitment |
plotMeasurementsVExpectedCost(measured_decisions) |
Plot optimizer output vs. classical optimum |
normalize(val) / unNormalize(val) |
Cost normalization helpers |
Solves the UC problem by keeping scenarios in superposition in the PDF register, measuring gas decisions only. This is the approach that makes QAE applicable.
from optimizer_utils import PowerSystem_1Bus
from dense_optimizer import Optimizer_Dense
system = PowerSystem_1Bus(...)
opt = Optimizer_Dense(system, encoding='binary')Annealing solve methods — all return a counts dict {bitstring: probability}:
| Method | Schedule | Notes |
|---|---|---|
solveAnnealing(time, method, num_meas, penalty) |
Linear, single-phase | QUBO or quadratic penalty; all-to-one transpile |
solveAnnealingAlternating(total_time, time_steps, ...) |
Alternating operator | QAOA-style; supports phase='PEN'/'COST' and mixer='X'/'SWAP'
|
solveAnnealingDenseStochHam(total_time, t2, t1_steps, t2_steps, ...) |
Two-phase stochastic | Phase 1: annealing; Phase 2: stochastic Hamiltonian |
solveAnnealingDenseStochHam2Layers(t1_time, t1_steps, t2_time, t2_steps, ...) |
Two-phase layered | Like above with explicit layer separation |
solveAnnealingStochHam(total_time, t_steps, ...) |
Stochastic Hamiltonian | Full stochastic Hamiltonian throughout |
solveECAnnealing(total_time, gamma, ...) |
EC (ensemble copy) | Multiple copies with weight |
solveThreePhaseAnnealing(total_time, ...) |
Three-phase | Annealing → PDF projection → re-annealing |
solveThreePhaseAnnealing(...) |
Three-phase | Phase 1: DQA; Phase 2: fixed; Phase 3: sampling with PDF reinit |
solveTwoPhaseAnnealing(total_time, time_2, ...) |
Two-phase sampling | Phase 1: build state; Phase 2: measure-and-reinit loop |
Circuit-building helpers:
| Method | Description |
|---|---|
priceOperator(amp) |
Gas cost phase operator |
scenarioOperator(amp) |
Scenario (wind) cost phase operator |
penaltyOperator(amp, penalty) |
Demand-constraint penalty |
initializePDF(inverse) |
Prepare (or un-prepare) the PDF register |
dickeState() |
Prepare Dicke state on decision qubits |
setVariables(vars, vals) |
Initialize a register to specific values |
Post-processing:
| Method | Description |
|---|---|
calcCostPerScenario(counts) |
Compute expected cost per scenario from measurements |
calcGasExpectationValue(counts) |
Expected value of gas commitment |
getGasCounts(counts) |
Marginalise over scenarios, return gas-only counts |
getDecision(gas_decisions) |
Extract optimal gas commitment from measurement histogram |
Alternative formulation: each scenario gets its own set of wind+slack qubits in the register (no superposition over scenarios). Simpler to understand, but exponentially larger register.
opt = Optimizer_Expanded(system, encoding='binary')Script for convergence experiments. Run directly:
python experiment_convergence.py --n 4 --m 3Generates QAE convergence data and saves to pickle for downstream analysis notebooks.
module load qiskit/aer-gpu
# Provides: Python 3.11, Qiskit 2.2.3, qiskit-aer-gpu 0.15.1 (CUDA 12.4, H100),
# qiskit-algorithms 0.4.0, matplotlib 3.10.8pip install qiskit>=2.0 qiskit-aer qiskit-algorithms matplotlib scipy numpy
# For resource estimation:
pip install pytket qualtranNote:
QuantumExpectedValueFunctionProject/code was written for an older Qiskit API. Before running it, apply these import fixes:# OLD NEW from qiskit import Aer → from qiskit_aer import Aer from qiskit.tools.visualization import … → from qiskit.visualization import … from qiskit.extensions import Initialize → from qiskit.circuit.library import Initialize
import numpy as np
from qiskit_impl.binary_optimizer import BinaryNestedOptimizer
# Problem: 2 wind turbines, total demand = 2 MW (both turbines needed).
# Each turbine produces exactly 1 MW when ON and wind is available (ξ_j = 1).
n = 2
pdf = {(0,0): 0.25, # no wind for either turbine
(0,1): 0.25, # wind only for turbine 1
(1,0): 0.25, # wind only for turbine 0
(1,1): 0.25} # wind available for both turbines
cy = np.linspace(0.01, 0.10, n) # wind dispatch costs per turbine ($/MWh)
cr = 1.0 # recourse cost for unmet demand ($/MWh)
cx = [0.4] # gas generator commitment cost ($/MW)
bno = BinaryNestedOptimizer(cx, cy, cr, pdf, demand=2, is_uniform=True)
# 1. Classical brute-force reference: φ(x) for ALL x in {0, 1, 2}
# ev_true[x] = E_ξ[min_y q(y,ξ)] when gas covers x MW, wind covers (demand-x) MW.
# Total cost for each x: 0.4*x + ev_true[x] → pick the x with minimum total cost.
# x=0: $0.00 + $1.075 = $1.075
# x=1: $0.40 + $0.300 = $0.700 ← optimal
# x=2: $0.80 + $0.000 = $0.800
ev_true = bno.brute_force_wind_demand_expectation_values()
print(f"Classical φ(x) for all x: {ev_true}")
# 2. DQA quantum estimate of φ(x=1) — the optimal first-stage decision
# 'wind_demand' encodes the first-stage gas decision x implicitly:
# wind_demand = demand - x
# wind_demand=2 → x=0 (gas off, wind covers all 2 MW)
# wind_demand=1 → x=1 (gas covers 1 MW, wind covers 1 MW) ← optimal
# wind_demand=0 → x=2 (gas covers everything, no wind needed)
# With wind_demand=1, Dicke state = (1/√2)(|10⟩ + |01⟩): turbine 0 or 1 on.
# DQA shifts amplitude toward turbine 0 (cheaper at $0.05 vs $0.10).
qc = bno.adiabatic_evolution_circuit(wind_demand=1, time=100, time_steps=4, norm=10)
counts = bno.execute_optimizer(qc, num_meas=8192)
exp_val = bno.process_expectation_value_optimizer(wind_demand=1, counts=counts)
print(f"DQA estimate of φ(x=1): {exp_val:.4f} (classical: {ev_true[1]:.4f})")
# 3. Full QAE pipeline (more accurate estimate of φ(x=1))
oracle = bno.single_oracle_sin_inconstraint(c=0.5, norm=10)
orc_inv = bno.single_oracle_sin_inconstraint(c=0.5, norm=10, inverse=True)
uopt = bno.adiabatic_evolution_circuit(wind_demand=1, time=100, time_steps=4, norm=10)
qae_qc = bno.canonical_qae(uopt, oracle, m=3, c=0.5, norm=10)
amplitude = bno.execute_qae(qae_qc, m=3, num_meas=8192)
print(f"QAE amplitude estimate: {amplitude:.4f}")import sys
sys.path.insert(0, 'QuantumExpectedValueFunctionProject')
from optimizer_utils import PowerSystem_1Bus
from dense_optimizer import Optimizer_Dense
system = PowerSystem_1Bus(
gas_costs=[0.4], wind_costs=[0.05, 0.10],
decision_levels=4, undersatisfied_cost=1.0,
demand=2,
pdf={(0,1):0.5, (1,0):0.5},
)
opt = Optimizer_Dense(system, encoding='binary')
counts = opt.solveAnnealingAlternating(
total_time=10, time_steps=5, penalty=1, num_meas=2000,
init_cond='XGS', mixer='X', phase='PEN',
)
decision = opt.getDecision(opt.getGasCounts(counts))
print("Gas commitment decision:", decision)| Notebook | Location | Purpose |
|---|---|---|
qae_example.ipynb |
qiskit_impl/ |
Step-by-step walkthrough of DQA + canonical QAE |
four_binary_turbines_uniform_distribution_tutorial.ipynb |
QuantumExpectedValueFunctionProject/ |
Tutorial: 4 turbines, uniform PDF |
ExpValFun_QAOA_QAE_computations.ipynb |
QuantumExpectedValueFunctionProject/ |
Full computation pipeline for QAOA+QAE |
ExpValFun_QAOA_QAE_figures.ipynb |
QuantumExpectedValueFunctionProject/ |
Publication-quality figures |
simple_UC_model.ipynb |
QuantumExpectedValueFunctionProject/ |
Simplest possible UC model |
SED_classical_solution.ipynb |
QuantumExpectedValueFunctionProject/ |
Classical stochastic ED reference |
On NREL HPC:
module load qiskit/aer-gpu
jupyter labBranch gpuOptim contains detailed comments marking every location that needs changing to run on the H100 GPU. No code is changed — only comments.
Replace every Aer.get_backend(...) call with:
from qiskit_aer import AerSimulator
simulator = AerSimulator(method='statevector', device='GPU')Marked with # GPU-SWAP TIER 1 in:
qiskit_impl/binary_optimizer.py—execute_optimizer(),execute_optimizer_oracle(),execute_qae()QuantumExpectedValueFunctionProject/dense_optimizer.py— all 8solveAnnealing*()methods
Inside DQA time-step loops, sub-circuits are currently rebuilt with numeric angles each iteration. The fix — demonstrated in bayesianQC/optimize_10epoch_performance.py:
# Build symbolic template ONCE
from qiskit.circuit import Parameter
s = Parameter('s')
cost_tmpl = cost_operator_symbolic(s, ...) # symbolic angles
# Transpile ONCE to GPU target
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
pm = generate_preset_pass_manager(optimization_level=1, backend=gpu_simulator)
cost_t = pm.run(cost_tmpl)
# Inside loop: only bind numbers (cheap)
for t in range(time_steps):
s_val = (dt * t) / time
qc.append(cost_t.assign_parameters({s: s_val}), qubits)Marked with # PARAM-TRANSPILE in binary_optimizer.py and dense_optimizer.py.
simulator = AerSimulator(
method='statevector', device='GPU',
cuStateVec_enable=True, # NVIDIA cuStateVec via cuquantum-cu12
blocking_enable=True, # auto-chunk if VRAM insufficient
blocking_qubits=23,
)from pytket import Circuit as TKCircuit
from qiskit_impl.resource_estimator import Resource_Estimator, qasm_to_clifford_and_t
from qiskit import QuantumCircuit
# Build the circuit you want to estimate
qc = QuantumCircuit(10)
# ... add gates ...
# Convert to pytket (for gate counting)
from qiskit.qasm2 import dumps
import pytket.extensions.qiskit as tkqk
tk_circ = tkqk.qiskit_to_tk(qc)
args = {
'num_algo_qubits' : qc.num_qubits,
'circuit' : tk_circ,
'num_meas' : 1000,
'error_budget' : 1e-3,
'architecture' : 'superconducting',
'optimistic' : False,
'R': 2, 'Q_0': 1e-3, 'P_r': 1e-6, 'F_copies': 1,
}
result = Resource_Estimator(args).bev_resource_estimator()
print(result)Based on the Beverland et al. (2022) surface-code resource model (arXiv:2211.07629).
- Each
file.pyhas a correspondingtest_file.py. Run tests withpython test_file.py. - In-development functions live as
_TEST_functionname()inside the source file, called frommain(). They are moved totest_file.pywhen the feature is complete. - Experiments go in Jupyter notebooks that import from the source modules — not in source modules themselves.
- Branch
nilesh/comments— adds explanatory comments and Qiskit 2.x compatibility fixes. - Branch
gpuOptim— adds GPU porting strategy comments (no code changes).
- Sævarsson et al., "Quantum algorithms for the two-stage stochastic unit commitment problem", arXiv:2402.15029 (2024)
- Brassard et al., "Quantum amplitude amplification and estimation", AMS Contemporary Mathematics 305 (2002)
- Beverland et al., "Assessing requirements to scale to practical quantum advantage", arXiv:2211.07629 (2022)
- Farhi et al., "A Quantum Approximate Optimization Algorithm", arXiv:1411.4028 (2014)
- Bartschi & Eidenbenz, "Deterministic Preparation of Dicke States", FCT 2019