Quantum computing is still at an early stage. The hardware has a long way to go before it can contribute to real-world problems, and the advantage it promises is over specific classes of task rather than computation in general. Even so, governments and companies continue to invest, preparing for the point at which the technology becomes usable. Chemistry, finance and energy are among the fields expected to benefit.
Combinatorial optimization is one of those classes of task, and it is the one we work on. Deciding which power plants to turn on to meet demand is a small example: the number of possible combinations grows exponentially with the number of plants, so even a system of 50 plants has more than a quadrillion candidate solutions to consider in the worst case. Problems of this shape are everywhere in power systems planning.
Getting such a model onto a quantum computer is not straightforward, for two separate reasons.
The first is the format. These machines accept problems in a very specific form: minimize a quadratic function over binary variables, with no constraints. That format is called QUBO (Quadratic Unconstrained Binary Optimization), and it is far more restrictive than the general-purpose modeling most practitioners use. It is not tied to a single vendor, either. Quantum annealers consume it directly. Gate-based quantum computers reach it through variational algorithms such as QAOA and VQE, which encode the problem as a Hamiltonian whose ground state is the answer. A family of physics-inspired classical devices works on the same principle: coherent Ising machines, simulated bifurcation machines, digital annealers. So do conventional heuristics like simulated annealing and parallel tempering, which need no special hardware at all.
The second is fragmentation. Because the field is young, each type of quantum machine has its own way of receiving problems and returning results. Classical computers converged on this long ago (different brands and operating systems still give you broadly the same experience), but nothing equivalent exists here yet. Someone who already knows how to write a QUBO would still have to learn a separate interface for every platform they wanted to try.
Each word in the name is a restriction, and each restriction has to be paid for somewhere.
Binary. Variables can only take the values 0 or 1. Anything else (integers, continuous quantities like the output of a power plant in MW) has to be rebuilt out of bits, which adds variables and, for continuous quantities, gives up precision.
Unconstrained. There is nowhere in the model to put a constraint. Each one is folded into the objective instead, as a penalty term that adds cost when it is violated and nothing when it is satisfied. How much cost is a decision in itself.
Quadratic. No term may involve more than two variables multiplied together. Penalty terms routinely come out with a higher degree than that, so they have to be broken back down by introducing auxiliary variables.
The final O is just optimization. Each of these steps has design choices that affect the quality of the model you end up with, and doing them by hand is slow and easy to get wrong.
QUBO.jl addresses both problems for models written in JuMP, Julia’s optimization modeling language. We develop it at PSR together with researchers from UFRJ, Purdue’s SECQUOIA group, and USRA at NASA Ames; the paper describing it appeared in Optimization Methods and Software in 2026.
What a QUBO is
You have a set of binary variables. Setting any one of them to 1 has a cost. Setting certain pairs to 1 together has an additional cost or bonus. Find the assignment that minimizes the total:
Q is a matrix of pairwise coefficients, with the diagonal holding the individual costs. There is nothing else in the model: no constraints, no other variable types.
Two facts make this restricted format worth caring about.
It is NP-hard when Q is not diagonal, so it is expressive enough to encode many combinatorial problems: scheduling, graph partitioning, routing, portfolio selection. Formulations for a large catalogue of NP problems are known.
And it is equivalent to the Ising model from statistical physics. Substituting s = 2x - 1 maps binary variables onto ±1 spins, turning the problem into finding the ground state of a system of interacting magnets. That is the reason the platforms listed above converge on this format: annealing hardware relaxes into low-energy states by construction, and variational algorithms on gate-based machines search for the ground state of the corresponding Hamiltonian. A QUBO is what an optimization problem looks like once it has been rewritten as physics.
The compiler framing
It helps to think of a QUBO instance as the assembly code of this hardware, and ToQUBO.jl, the reformulation package in the ecosystem, as the compiler that produces it. A JuMP model goes in, a QUBO comes out, and we keep the mapping between the two in memory so results can be reported in terms of the original variables.
ToQUBO.jl is implemented as a MathOptInterface layer, the same mechanism used by JuMP extensions like Dualization.jl and QuadraticToBinary.jl. The QUBO it produces is itself a valid MOI model, so it can go to a quantum annealer or to a conventional non-convex MIQP solver like Gurobi without changing anything upstream.
A worked example: thermal dispatch
The clearest way to show what compilation involves is to do it by hand on the problem from the introduction.
Take a simplified thermal dispatch: N plants, each with an operating cost cᵢ in $/MW and a maximum output Gᵢᵐᵃˣ. Meet a demand D as cheaply as possible.
With three plants (costs of 1, 2 and 3 $/MW, capacities of 2, 3 and 5 MW, and a demand of 6 MW), this is:
You can solve it in your head. Dispatch the cheap plants first: g = (2, 3, 1), for a cost of 11. That makes it a good test case, because we know what the answer should be.
In JuMP, with QUBO.jl as the solver:
using JuMP
using QUBO
c = [1, 2, 3] # operating cost, $/MW
G = [2, 3, 5] # capacity, MW
D = 6 # demand, MW
model = Model(() -> ToQUBO.Optimizer(ExactSampler.Optimizer))
@variable(model, 0 <= g[i in 1:3] <= G[i], Int)
@objective(model, Min, c' * g)
@constraint(model, sum(g) == D)
optimize!(model)Integer variables and an equality constraint: two of the three restrictions above are violated on sight. Here is what it takes to remove them.
Stage 1: penalization
The constraint moves into the objective as a squared violation term, weighted by a penalty factor ρ:
The term is zero when supply matches demand and positive otherwise. How large does ρ need to be? Consider shutting off the last MW of the expensive plant: that saves 3 in operating cost and adds ρ in penalty. So ρ has to exceed 3, the largest cost coefficient in the problem. Otherwise the compiled model prefers a blackout to a correct answer. Too large and the cost terms get swamped by penalty terms, which hurts conditioning. Choosing ρ per constraint is one of the things ToQUBO.jl does automatically, and one you can override with a hint when you know better.
Stage 2: variable encoding
Each integer variable now has to be rebuilt from bits. For an integer I ∈ [a, b] we use n = ⌈log₂(b − a) + 1⌉ binary variables:
The leading term corrects the top bit so the range comes out exactly right instead of overshooting to the next power of two. Applied to the three plants:
| Variable | Range | Bits | Expansion |
|---|---|---|---|
g₁ | [0, 2] | 2 | g₁₁ + g₁₂ |
g₂ | [0, 3] | 3 | g₂₁ + 2g₂₂ |
g₃ | [0, 5] | 4 | g₃₁ + 2g₃₂ + 4g₃₃ − 2g₃₄ |
Note g₃₄ carrying a negative weight, and note that g₂₃ drops out entirely with a coefficient of zero: the formula allocated a bit that this particular range doesn’t need.
Substituting everything back:
Stage 3: quadratization
Expanding that square produces only pairwise products, so this particular model is already of degree two and no further work is needed.
That is not true in general. Penalty functions for other constraint families (logical constraints, quadratic constraints) routinely produce terms of degree three or higher, which have to be reduced by introducing auxiliary variables. We ship two standard single-term reduction methods, and the interface is extensible through Julia’s multiple dispatch if you want to add your own.
What came out
Three integer variables and one constraint went in. What comes out is a dense 9×9 Q matrix over nine binary variables, with entries ranging from ρ to 32ρ, one variable that does nothing, and off-diagonal coupling between nearly every pair.
This is the real cost of the format. A problem you can solve by inspection becomes a dense quadratic nobody would want to write out, and that is the argument for having a compiler generate it rather than a person.
Encoding integers, and why coefficient size matters
The dispatch example used binary place-value encoding, with weights 1, 2, 4, 8. It is the most economical choice in variable count: an integer up to n needs only O(log n) bits. But the largest weight grows as O(n), and in the compiled model above that already produced a 32× spread between the smallest and largest entries of Q from a problem with three plants.
Coefficient magnitude is not a cosmetic concern. On annealing hardware, coefficients correspond to physical quantities (magnetic field strengths, coupling terms), and the available precision is finite. When the ratio between the largest and smallest coefficients is wide, the small ones fall below the noise floor and the machine effectively solves a different problem than the one you specified. Scale the dispatch example up to a real system with hundreds of plants and the spread gets much worse.
The usual alternative is unary encoding: n variables all of weight 1. Coefficients stay at O(1), but the variable count is linear, and on hardware with a few thousand qubits that ceiling arrives quickly. Other schemes (one-hot, domain-wall, bounded-coefficient) sit at different points on the same trade-off.
| Encoding | Variables | Largest coefficient |
|---|---|---|
| Binary | O(log n) | O(n) |
| Unary | O(n) | O(1) |
| One-hot | O(n) | O(n) |
| Domain-wall | O(n) | O(n) |
| Bounded-coefficient | O(n) | O(1) |
There is no globally correct choice here; it depends on the size of the variable’s range and on how much precision the target device actually has. ToQUBO.jl implements all of these and lets you set the encoding per variable, so the same JuMP model can be recompiled for different hardware without being rewritten.
The rest of the ecosystem
QUBODrivers.jl is the answer to the second problem from the introduction. It defines a common MathOptInterface-compliant API for QUBO samplers, so that a piece of annealing hardware appears to JuMP as an ordinary solver and switching between platforms means changing one line. We have published wrappers for D-Wave’s quantum and simulated annealers, IBM hardware via Qiskit (QAOA and VQE), Los Alamos’ QuantumAnnealing.jl simulator, coherent Ising machine simulators, NASA’s PySA parallel tempering solver, and the MQLib heuristics library. New wrappers are written against a macro-based setup interface rather than against MOI internals, which keeps them short: if you build annealing hardware, connecting it is roughly a page of Julia.
QUBOTools.jl handles I/O and analysis: conversion between the file formats different vendors use (bqpjson, D-Wave’s QUBO and Qubist formats, others), reading public benchmark instance databases, computing metrics like time-to-solution and success rate, and plotting recipes for model density and sample distributions.
Benchmarks
We compared model-building time against PyQUBO, qubovert, Qiskit, and OpenQAOA, following PyQUBO’s own benchmarking setup: traveling salesman instances from 25 to 10,000 variables, and number partitioning from 5 to 1,000. ToQUBO.jl builds models faster than PyQUBO and substantially faster than the other three. The benchmark code is public, so the numbers can be reproduced.
On features, ToQUBO.jl covers every encoding and constraint type the others support, and is the only tool in the comparison whose encoding routines also apply to continuous variables rather than integers alone. The coverage table in the paper has the full breakdown.
Part of the performance difference comes down to the language. PyQUBO’s reformulation core had to be rewritten from Python to C++ for speed, which is exactly the two-language problem Julia is designed to avoid. That was a large part of why we chose Julia in the first place.
Where this stands
The caveat from the opening still applies at the end: none of this makes quantum hardware competitive today, and we are not claiming otherwise.
What the tooling does provide is useful regardless. QUBO is the input format for GPU and FPGA-based annealers and for mature classical heuristics like MQLib, none of which are quantum. Compiling a model once makes all of them reachable, and gives a way to run the same instance across them, which is a prerequisite for judging whether any of it works. That comparison is the part we are most interested in, and it is where the next round of work on the ecosystem is going.
- Repository: github.com/JuliaQUBO/QUBO.jl
- Documentation: juliaqubo.github.io/QUBO.jl
- Paper: QUBO.jl: A Julia Ecosystem for Quadratic Unconstrained Binary Optimization, Optimization Methods and Software, 2026. DOI · arXiv preprint
