Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

QUBO and Ising Models (Julia)

QUBO & Ising ModelsΒΆ

David E. Bernal Neira
Davidson School of Chemical Engineering, Purdue University

Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University
Computer Science & Systems Engineering Program, Federal University of Rio de Janeiro
PSR Energy Consulting & Analytics

JoΓ£o Victor Souza
Computer Engineering Department, Military Institute of Engineering (IME)

Open In Colab SECQUOIA

Environment setup

For local execution from the repository root, run uv sync --group qubo and make setup-julia NOTEBOOK=notebooks_jl/2-QUBO.ipynb before launching Jupyter. This notebook reuses the repo-local Python environment for the D-Wave Ocean stack instead of relying on Julia’s CondaPkg resolver.

This notebook should open directly with the Julia runtime in Google Colab. If it does not, switch the runtime before running the next cell. The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/envs/2-QUBO/Project.toml, install the Python D-Wave Ocean packages into the Colab Python runtime, and print progress for clone, activate, instantiate, and package-load steps. Package precompilation still happens as needed during setup. The first Colab run can still take several minutes.

function load_notebook_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QuIP", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QuIP", "scripts", "notebook_bootstrap.jl"),
    )

    for candidate in candidates
        if isfile(candidate)
            include(candidate)
            return nothing
        end
    end

    in_colab = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP") || isdir(joinpath("/content", "sample_data"))
    if in_colab
        repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning SECQUOIA/QuIP into $repo_dir")
            run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
        end
        include(joinpath(repo_dir, "scripts", "notebook_bootstrap.jl"))
        return nothing
    end

    error("Could not locate scripts/notebook_bootstrap.jl from $(pwd()).")
end

load_notebook_bootstrap()

BOOTSTRAP = QuIPNotebookBootstrap.bootstrap_notebook("2-QUBO")

QUIP_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab

if !@isdefined(Karnak)
    using Karnak
end

if !@isdefined(LinearAlgebra)
    using LinearAlgebra
end

if !@isdefined(Graphs)
    using Graphs
end

if !@isdefined(JuMP)
    using JuMP
end

if !@isdefined(QUBO)
    using QUBO
end

if !@isdefined(Plots)
    using Plots
end

if !@isdefined(HiGHS)
    using HiGHS
end

if !@isdefined(DWave)
    using DWave
end

if !@isdefined(Luxor)
    using Luxor
end

QUBO and Ising Models (Julia)ΒΆ

This notebook explains the basics of QUBO modeling. In order to implement the different QUBO formulations we will use JuMP, and then solve them using neal’s implementation of simulated annealing. We will also leverage QUBO.jl to translate constraint satisfaction problems to QUBOs. Finally, we will use Graphs.jl for network models and graph visualizations.

Problem statementΒΆ

We define a QUBO as the following optimization problem:

min⁑x∈{0,1}nβˆ‘(ij)∈E(G)Qijxixj+βˆ‘i∈V(G)Qiixi+Ξ²=min⁑x∈{0,1}nxβ€²Qx+Ξ²\min_{x \in \{0,1 \}^n} \sum_{(ij) \in E(G)} Q_{ij}x_i x_j + \sum_{i \in V(G)}Q_{ii}x_i + \beta = \min_{x \in \{0,1 \}^n} \mathbf{x}' \mathbf{Q} \mathbf{x} + \beta
(1)

where we optimize over binary variables x∈{0,1}nx \in \{ 0,1 \}^n, on a constrained graph G(V,E)G(V,E) defined by a weighted adjacency matrix Q\mathbf{Q}. We also include an arbitrary offset β\beta.

ExampleΒΆ

Suppose we want to solve the following problem via QUBO

min⁑x2x0+4x1+4x2+4x3+4x4+4x5+5x6+4x7+5x8+6x9+5x10s.t.[100111011110101011011100101011111]x=[111] x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} & 2x_0+4x_1+4x_2+4x_3+4x_4+4x_5+5x_6+4x_7+5x_8+6x_9+5x_{10} \\ \textrm{s.t.} & \begin{bmatrix} 1 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1\\ 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1\\ 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}\mathbf{x}= \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix} \\ ~ & \mathbf{x} \in \{0,1 \}^{11} \end{array}
(2)

First we would write this problem as an unconstrained one by penalizing the linear constraints as quadratics in the objective. Let’s first define the problem parameters.

A = [
    1 0 0 1 1 1 0 1 1 1 1
    0 1 0 1 0 1 1 0 1 1 1
    0 0 1 0 1 0 1 1 1 1 1
]
b = [1, 1, 1]
c = [2, 4, 4, 4, 4, 4, 5, 4, 5, 6, 5];

In order to define the Q\mathbf{Q} matrix, we first write the problem

min⁑xcβ€²xs.t.Ax=bΒ x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} &\mathbf{c}' \mathbf{x} \\ \textrm{s.t.} & \mathbf{A}\mathbf{x} = \mathbf{b} \\ ~ & \mathbf{x} \in \{0,1 \}^{11} \end{array}
(3)

as follows:

min⁑xcβ€²x+ρ(Axβˆ’b)β€²(Axβˆ’b)s.t.x∈{0,1}11\begin{array}{rl} \displaystyle% \min_{\mathbf{x}} & \mathbf{c}' \mathbf{x} + \rho (\mathbf{A}\mathbf{x}-\mathbf{b})' (\mathbf{A}\mathbf{x}-\mathbf{b}) \\ \textrm{s.t.} & \mathbf{x} \in \{0,1 \}^{11} \end{array}
(4)

Exploiting the fact that x2=xx^2=x for x∈{0,1}x \in \{0,1\}, we can make the linear terms appear in the diagonal of the Q\mathbf{Q} matrix.

ρ(Axβˆ’b)β€²(Axβˆ’b)=ρ(xβ€²(Aβ€²A)xβˆ’2(Aβ€²b)x+bβ€²b)\rho(\mathbf{A}\mathbf{x}-\mathbf{b})'(\mathbf{A}\mathbf{x}-\mathbf{b}) = \rho( \mathbf{x}'(\mathbf{A}'\mathbf{A}) \mathbf{x} - 2(\mathbf{A}'\mathbf{b}) \mathbf{x} + \mathbf{b}'\mathbf{b} )
(5)

For this problem in particular, one can prove that a reasonable penalization factor is given by ρ=βˆ‘i=1n∣ci∣+Ο΅\rho = \sum_{i=1}^n |c_i| + \epsilon with Ο΅>0\epsilon > 0.

Ο΅ = 1
ρ = sum(abs, c) + ϡ
Q = diagm(c) + ρ * (A'A - 2 * diagm(A'b))
β = ρ * b'b

display(Q)
println(Ξ²)
11Γ—11 Matrix{Int64}: -46 0 0 48 48 48 0 48 48 48 48 0 -44 0 48 0 48 48 0 48 48 48 0 0 -44 0 48 0 48 48 48 48 48 48 48 0 -92 48 96 48 48 96 96 96 48 0 48 48 -92 48 48 96 96 96 96 48 48 0 96 48 -92 48 48 96 96 96 0 48 48 48 48 48 -91 48 96 96 96 48 0 48 48 96 48 48 -92 96 96 96 48 48 48 96 96 96 96 96 -139 144 144 48 48 48 96 96 96 96 96 144 -138 144 48 48 48 96 96 96 96 96 144 144 -139
144

We can visualize the graph that defines this instance using the Q\mathbf{Q} matrix as the adjacency matrix of a graph.

G = SimpleGraph(Q)
println(G)
SimpleGraph{
Int64}(58, 
[[1, 4, 5, 6, 8, 9, 10, 11], [2, 4, 6, 7, 9, 10, 11], [3, 5, 7, 8, 9, 10, 11], [1, 2, 4, 5, 6, 7, 8, 9, 10, 11], [1, 3, 4, 5, 6, 7, 8, 9, 10, 11], [1, 2, 4, 5, 6, 7, 8, 9, 10, 11], [2, 3, 4, 5, 6, 7, 8, 9, 10, 11], [1, 3, 4, 5, 6, 7, 8, 9, 10, 11], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]])

@drawsvg(
    begin
        sethue("black")
        background("white")
        drawgraph(
            G;
            margin       = 80,
            vertexlabels = 1:11,
        )
    end,
    400, 400
)
Loading...

Let’s define a QUBO model and then solve it using enumeration and D-Wave’s simulated annealing (and, eventually, quantum annealing too!).

# Define empty model
qubo_model = Model()

# Define the variables
@variable(qubo_model, x[1:11], Bin)

# Define the objective function
@objective(qubo_model, Min, x' * Q * x + Ξ²)

# Print the model
print(qubo_model)
Loading...

Since the problem is relatively small (11 variables, 211=20482^{11}=2048 combinations), we can afford to enumerate all the solutions.

# Here we solve the optimization problem by exact enumeration.
set_optimizer(qubo_model, ExactSampler.Optimizer)

optimize!(qubo_model)

qubo_x = round.(Int, value.(x))
qubo_objective = objective_value(qubo_model)

# Display solution of the problem
println(solution_summary(qubo_model))
println("* objective value = $(round(qubo_objective; digits=1))")
println("* x = $qubo_x")
solution_summary(; result = 1, verbose = false)
β”œ solver_name          : Exact Sampler
β”œ Termination
β”‚ β”œ termination_status : LOCALLY_SOLVED
β”‚ β”œ result_count       : 2048
β”‚ β”” raw_status         : optimal
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : 
FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”œ objective_value      : 
5.00000e+00
β”‚ β”” dual_objective_value : 1.44000e+02
β”” Work counters
  β”” solve_time (sec)   : 1.65687e-01
* objective value = 5.0
* x = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]

plot(QUBOTools.EnergyFrequencyPlot(unsafe_backend(qubo_model)))
Loading...

Let’s now solve this QUBO via traditional Integer Programming.

qubo_ilp_model = Model()

@variable(qubo_ilp_model, x[1:11], Bin)
@variable(qubo_ilp_model, y[1:11, 1:11], Bin)

@objective(
    qubo_ilp_model,
    Min,
    sum(Q[i,j] * (i == j ? x[i] : y[i,j]) for i=1:11, j=1:11) +  Ξ²
)

@constraint(qubo_ilp_model, c1[i=1:11,j=1:11;i!=j], y[i,j] >= x[i] + x[j] - 1)
@constraint(qubo_ilp_model, c2[i=1:11,j=1:11;i!=j], y[i,j] <= x[i])
@constraint(qubo_ilp_model, c3[i=1:11,j=1:11;i!=j], y[i,j] <= x[j])

println(qubo_ilp_model)
Min -46 x[1] + 48 y[1,4] + 48 y[1,5] + 48 y[1,6] + 48 y[1,8] + 48 y[1,9] + 48 y[1,10] + 48 y[1,11] - 44 x[2] + 48 y[2,4] + 48 y[2,6] + 48 y[2,7] + 48 y[2,9] + 48 y[2,10] + 48 y[2,11] - 44 x[3] + 48 y[3,5] + 48 y[3,7] + 48 y[3,8] + 48 y[3,9] + 48 y[3,10] + 48 y[3,11] + 48 y[4,1] + 48 y[4,2] - 92 x[4] + 48 y[4,5] + 96 y[4,6] + 48 y[4,7] + 48 y[4,8] + 96 y[4,9] + [[...45 terms omitted...]] + 96 y[9,4] + 96 y[9,5] + 96 y[9,6] + 96 y[9,7] + 96 y[9,8] - 139 x[9] + 144 y[9,10] + 144 y[9,11] + 48 y[10,1] + 48 y[10,2] + 48 y[10,3] + 96 y[10,4] + 96 y[10,5] + 96 y[10,6] + 96 y[10,7] + 96 y[10,8] + 144 y[10,9] - 138 x[10] + 144 y[10,11] + 48 y[11,1] + 48 y[11,2] + 48 y[11,3] + 96 y[11,4] + 96 y[11,5] + 96 y[11,6] + 96 y[11,7] + 96 y[11,8] + 144 y[11,9] + 144 y[11,10] - 139 x[11] + 144

Subject to
 
c1[1,2] : -x[1] - x[2] + y[1,2] β‰₯ -1
 c1[1,3] : -x[1] - x[3] + y[1,3] β‰₯ -1
 c1[1,4] : -x[1] - x[4] + y[1,4] β‰₯ -1
 c1[1,5] : -x[1] - x[5] + y[1,5] β‰₯ -1
 c1[1,6] : -x[1] - x[6] + y[1,6] β‰₯ -1
 c1[1,7] : -x[1] - x[7] + y[1,7] β‰₯ -1
 c1[1,8] : -x[1] - x[8] + y[1,8] β‰₯ -1
 c1[1,9] : -x[1] - x[9] + y[1,9] β‰₯ -1
 c1[1,10] : -x[1] - x[10] + y[1,10] β‰₯ -1
 c1[1,11] : -x[1] - x[11] + y[1,11] β‰₯ -1
 c1[2,1] : -x[1] - x[2] + y[2,1] β‰₯ -1
 c1[2,3] : -x[2] - x[3] + y[2,3] β‰₯ -1
 c1[2,4] : -x[2] - x[4] + y[2,4] β‰₯ -1
 c1[2,5] : -x[2] - x[5] + y[2,5] β‰₯ -1
 c1[2,6] : -x[2] - x[6] + y[2,6] β‰₯ -1
 c1[2,7] : -x[2] - x[7] + y[2,7] β‰₯ -1
 c1[2,8] : -x[2] - x[8] + y[2,8] β‰₯ -1
 c1[2,9] : -x[2] - x[9] + y[2,9] β‰₯ -1
 c1[2,10] : -x[2] - x[10] + y[2,10] β‰₯ -1
 c1[2,11] : -x[2] - x[11] + y[2,11] β‰₯ -1
 c1[3,1] : -x[1] - x[3] + y[3,1] β‰₯ -1
 c1[3,2] : -x[2] - x[3] + y[3,2] β‰₯ -1
 c1[3,4] : -x[3] - x[4] + y[3,4] β‰₯ -1
 c1[3,5] : -x[3] - x[5] + y[3,5] β‰₯ -1
 c1[3,6] : -x[3] - x[6] + y[3,6] β‰₯ -1
 c1[3,7] : -x[3] - x[7] + y[3,7] β‰₯ -1
 c1[3,8] : -x[3] - x[8] + y[3,8] β‰₯ -1
 c1[3,9] : -x[3] - x[9] + y[3,9] β‰₯ -1
 c1[3,10] : -x[3] - x[10] + y[3,10] β‰₯ -1
 c1[3,11] : -x[3] - x[11] + y[3,11] β‰₯ -1
 c1[4,1] : -x[1] - x[4] + y[4,1] β‰₯ -1
 c1[4,2] : -x[2] - x[4] + y[4,2] β‰₯ -1
 c1[4,3] : -x[3] - x[4] + y[4,3] β‰₯ -1
 c1[4,5] : -x[4] - x[5] + y[4,5] β‰₯ -1
 c1[4,6] : -x[4] - x[6] + y[4,6] β‰₯ -1
 c1[4,7] : -x[4] - x[7] + y[4,7] β‰₯ -1
 c1[4,8] : -x[4] - x[8] + y[4,8] β‰₯ -1
 c1[4,9] : -x[4] - x[9] + y[4,9] β‰₯ -1
 c1[4,10] : -x[4] - x[10] + y[4,10] β‰₯ -1
 c1[4,11] : -x[4] - x[11] + y[4,11] β‰₯ -1
 c1[5,1] : -x[1] - x[5] + y[5,1] β‰₯ -1
 c1[5,2] : -x[2] - x[5] + y[5,2] β‰₯ -1
 c1[5,3] : -x[3] - x[5] + y[5,3] β‰₯ -1
 c1[5,4] : -x[4] - x[5] + y[5,4] β‰₯ -1
 c1[5,6] : -x[5] - x[6] + y[5,6] β‰₯ -1
 c1[5,7] : -x[5] - x[7] + y[5,7] β‰₯ -1
 c1[5,8] : -x[5] - x[8] + y[5,8] β‰₯ -1
 c1[5,9] : -x[5] - x[9] + y[5,9] β‰₯ -1
 c1[5,10] : -x[5] - x[10] + y[5,10] β‰₯ -1
 c1[5,11] : -x[5] - x[11] + y[5,11] β‰₯ -1
[[...362 constraints skipped...]]
 y[6,7] binary
 y[7,7] binary
 y[8,7] binary
 y[9,7] binary
 y[10,7] binary
 y[11,7] binary
 y[1,8] binary
 y[2,8] binary
 y[3,8] binary
 y[4,8] binary
 y[5,8] binary
 y[6,8] binary
 y[7,8] binary
 y[8,8] binary
 y[9,8] binary
 y[10,8] binary
 y[11,8] binary
 y[1,9] binary
 y[2,9] binary
 y[3,9] binary
 y[4,9] binary
 y[5,9] binary
 y[6,9] binary
 y[7,9] binary
 y[8,9] binary
 y[9,9] binary
 y[10,9] binary
 y[11,9] binary
 y[1,10] binary
 y[2,10] binary
 y[3,10] binary
 y[4,10] binary
 y[5,10] binary
 y[6,10] binary
 y[7,10] binary
 y[8,10] binary
 y[9,10] binary
 y[10,10] binary
 y[11,10] binary
 y[1,11] binary
 y[2,11] binary
 y[3,11] binary
 y[4,11] binary
 y[5,11] binary
 y[6,11] binary
 y[7,11] binary
 y[8,11] binary
 y[9,11] binary
 y[10,11] binary
 y[11,11] binary


set_optimizer(qubo_ilp_model, HiGHS.Optimizer)

optimize!(qubo_ilp_model)

qubo_ilp_x = round.(Int, value.(x))
qubo_ilp_objective = objective_value(qubo_ilp_model)

println(solution_summary(qubo_ilp_model))
println("* objective value = $(round(qubo_ilp_objective; digits=1))")
println("* x = $qubo_ilp_x")
Running HiGHS 1.13.1 (git hash: 1d267d97c): Copyright (c) 2026 under Apache 2.0 license terms
Using BLAS: blastrampoline 
MIP has 330 rows; 132 cols; 770 nonzeros; 132 integer variables (132 binary)
Coefficient ranges:
  Matrix  [1e+00, 1e+00]
  Cost    [4e+01, 1e+02]
  Bound   [1e+00, 1e+00]
  RHS     [1e+00, 1e+00]
Presolving model
330 rows, 121 cols, 770 nonzeros  0s
165 rows, 66 cols, 385 nonzeros  0s
165 rows, 66 cols, 385 nonzeros  0s
Presolve reductions: rows 165(-165); columns 66(-66); nonzeros 385(-385) 
Objective function is integral with scale 1

Solving MIP model with:
   165 rows
   66 cols (66 binary, 0 integer, 0 implied int., 0 continuous, 0 domain fixed)
   385 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
     I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
     S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
     Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

 z       0       0         0   0.00%   -inf            144                Large        0      0      0         0     0.0s
         0       0         0   0.00%   -360.5          144              350.35%        0      0      5        13     0.0s
 L       0       0         0   0.00%   -140.9411765    5               2918.82%      228     18     10        73     0.1s

22.7% inactive integer columns, restarting
Model after restart has 135 rows, 51 cols (51 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 310 nonzeros

         0       0         0   0.00%   -124.8571429    5               2597.14%        6      0      0       118     0.1s
         0       0         0   0.00%   -124.8571429    5               2597.14%        6      6      2       131     0.1s

35.3% inactive integer columns, restarting
Model after restart has 99 rows, 33 cols (33 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 220 nonzeros

         0       0         0   0.00%   -86.5           5               1830.00%        8      0      0       157     0.1s
         0       0         0   0.00%   -86.5           5               1830.00%        8      8      5       170     0.1s

27.3% inactive integer columns, restarting
Model after restart has 60 rows, 24 cols (24 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 136 nonzeros

         0       0         0   0.00%   -77.16666667    5               1643.33%        5      0      0       183     0.1s
         0       0         0   0.00%   -77.16666667    5               1643.33%        5      5      3       199     0.2s
         1       0         1 100.00%   5               5                  0.00%       75      7      3       214     0.2s

Solving report
  Status            Optimal
  Primal bound      5
  Dual bound        5
  Gap               0% (tolerance: 0.01%)
  P-D integral      1.59638819802
  Solution status   feasible
                    5 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.15
  Max sub-MIP depth 2
  Nodes             1
  Repair LPs        0
  LP iterations     214
                    0 (strong br.)
                    125 (separation)
                    34 (heuristics)
solution_summary(; result = 1, verbose = false)

β”œ solver_name          : HiGHS
β”œ Termination
β”‚ β”œ termination_status : OPTIMAL
β”‚ β”œ result_count       : 1
β”‚ β”œ raw_status         : kHighsModelStatusOptimal
β”‚ β”” objective_bound    : 5.00000e+00
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”œ objective_value      : 5.00000e+00
β”‚ β”œ dual_objective_value : NaN
β”‚ β”” relative_gap         : 0.00000e+00
β”” Work counters
  β”œ solve_time (sec)   : 1.51395e-01
  β”œ simplex_iterations : 214
  β”œ barrier_iterations : -1
  β”” node_count         : 1
* objective value = 5.0
* x = [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]

Both the exact sampler and the HiGHS reformulation attain the minimum QUBO objective value of 5. In 1-based indexing, the two minimizing binary solutions have either x9=1x_9 = 1 or x11=1x_{11} = 1, with all other entries equal to 0.

Ising modelΒΆ

This notebook explains the basics of the Ising model. In order to implement the different Ising models we will use JuMP and D-Wave’s neal, for defining the Ising model and solving it with simulated annealing, respectively. When posing the problem as an integer program, we convert it to an equivalent binary QUBO formulation and then solve that MILP with the open-source HiGHS solver.

Problem statementΒΆ

We pose the Ising problem as the following optimization problem:

min⁑s∈{Β±1}nH(s)=min⁑s∈{Β±1}nβˆ‘(i,j)∈E(G)Ji,jsisj+βˆ‘i∈V(G)hisi+Ξ²\min_{s \in \{ \pm 1 \}^n} H(s) = \min_{s \in \{ \pm 1 \}^n} \sum_{(i, j) \in E(G)} J_{i,j}s_is_j + \sum_{i \in V(G)} h_is_i + \beta
(6)

where we optimize over spins s∈{±1}ns \in \{ \pm 1 \}^n, on a constrained graph G(V,E)G(V,E), where the quadratic coefficients are Ji,jJ_{i,j} and the linear coefficients are hih_i. We also include an arbitrary offset of the Ising model β\beta.

ExampleΒΆ

Suppose we have an Ising model defined from

h=[145.0122.0122.0266.0266.0266.0242.5266.0386.5387.0386.5],J=[00024242424242424240002402424242424240000240242424242400002448242448484800000242448484848000000242448484800000002448484800000000484848000000000727200000000007200000000000]Β andΒ Ξ²=1319.5h = \begin{bmatrix} 145.0 \\ 122.0 \\ 122.0 \\ 266.0 \\ 266.0 \\ 266.0 \\ 242.5 \\ 266.0 \\ 386.5 \\ 387.0 \\ 386.5 \end{bmatrix}, J = \begin{bmatrix} 0 & 0 & 0 & 24 & 24 & 24 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 24 & 0 & 24 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 0 & 24 & 0 & 24 & 24 & 24 & 24 & 24\\ 0 & 0 & 0 & 0 & 24 & 48 & 24 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 24 & 24 & 48 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 24 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 24 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 48 & 48 & 48\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 72 & 72\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 72\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ \end{bmatrix} \text{ and } \beta = 1319.5
(7)

Let’s solve this problem

h = [
    145.0,
    122.0,
    122.0,
    266.0,
    266.0,
    266.0,
    242.5,
    266.0,
    386.5,
    387.0,
    386.5,
]

J = [
    0  0  0 24 24 24  0 24 24 24 24
    0  0  0 24  0 24 24  0 24 24 24
    0  0  0  0 24  0 24 24 24 24 24
    0  0  0  0 24 48 24 24 48 48 48
    0  0  0  0  0 24 24 48 48 48 48
    0  0  0  0  0  0 24 24 48 48 48
    0  0  0  0  0  0  0 24 48 48 48
    0  0  0  0  0  0  0  0 48 48 48
    0  0  0  0  0  0  0  0  0 72 72
    0  0  0  0  0  0  0  0  0  0 72
    0  0  0  0  0  0  0  0  0  0  0
]

Ξ² = 1319.5

ising_model = Model()

@variable(ising_model, s[1:11], Spin)

@objective(ising_model, Min, s' * J * s + h' * s + Ξ²)

println(ising_model)
Min 24 s[4]*s[1] + 24 s[4]*s[2] + 24 s[5]*s[1] + 24 s[5]*s[3] + 24 s[5]*s[4] + 24 s[6]*s[1] + 24 s[6]*s[2] + 48 s[6]*s[4] + 24 s[6]*s[5] + 24 s[7]*s[2] + 24 s[7]*s[3] + 24 s[7]*s[4] + 24 s[7]*s[5] + 24 s[7]*s[6] + 24 s[8]*s[1] + 24 s[8]*s[3] + 24 s[8]*s[4] + 48 s[8]*s[5] + 24 s[8]*s[6] + 24 s[8]*s[7] + 24 s[9]*s[1] + 24 s[9]*s[2] + 24 s[9]*s[3] + 48 s[9]*s[4] + 48 s[9]*s[5] + 48 s[9]*s[6] + 48 s[9]*s[7] + 48 s[9]*s[8] + 24 s[10]*s[1] + 24 s[10]*s[2] + 24 s[10]*s[3] + 48 s[10]*s[4] + 48 s[10]*s[5] + 48 s[10]*s[6] + 48 s[10]*s[7] + 48 s[10]*s[8] + 72 s[10]*s[9] + 24 s[11]*s[1] + 24 s[11]*s[2] + 24 s[11]*s[3] + 48 s[11]*s[4] + 48 s[11]*s[5] + 48 s[11]*s[6] + 48 s[11]*s[7] + 48 s[11]*s[8] + 72 s[11]*s[9] + 72 s[11]*s[10] + 145 s[1] + 122 s[2] + 122 s[3] + 266 s[4] + 266 s[5] + 266 s[6] + 242.5 s[7] + 266 s[8] + 386.5 s[9] + 387 s[10] + 386.5 s[11] + 1319.5

Subject to
 
s[1] spin
 s[2] spin
 s[3] spin
 s[4] spin
 s[5] spin
 s[6] spin
 s[7] spin
 s[8] spin
 s[9] spin
 s[10] spin
 s[11] spin


set_optimizer(ising_model, ExactSampler.Optimizer)

optimize!(ising_model)

ising_s = round.(Int, value.(s))
ising_objective = objective_value(ising_model)

# Display solution of the problem
println(solution_summary(ising_model))
println("* objective value = $(round(ising_objective; digits=1))")
println("* s = $ising_s")
solution_summary(; result = 1, verbose = false)
β”œ solver_name          : Exact Sampler
β”œ Termination
β”‚ β”œ termination_status : LOCALLY_SOLVED
β”‚ β”œ result_count       : 2048
β”‚ β”” raw_status         : optimal
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”œ objective_value      : 5.00000e+00
β”‚ β”” dual_objective_value : 1.31950e+03
β”” Work counters
  β”” solve_time (sec)   : 7.16278e-04
* objective value = 5.0
* s = [-1, -1, -1, -1, -1, -1, -1, -1, 1, -1, -1]

plot(QUBOTools.EnergyFrequencyPlot(unsafe_backend(ising_model)))
Loading...
qubo_form = QUBOTools.qubo(unsafe_backend(ising_model), :dense; sense = :min)

n = QUBOTools.dimension(qubo_form)
Ξ± = QUBOTools.scale(qubo_form)
Ξ² = QUBOTools.offset(qubo_form)
Q = zeros(Float64, n, n)

for term in QUBOTools.linear_terms(qubo_form)
    i, value = term
    Q[i, i] = value
end

for term in QUBOTools.quadratic_terms(qubo_form)
    (i, j), value = term
    Q[i, j] = value
end
144.0
ising_ilp_model = Model()

@variable(ising_ilp_model, x[1:11], Bin)
@variable(ising_ilp_model, y[1:11, 1:11], Bin)

@objective(
    ising_ilp_model,
    Min,
    Ξ± * sum(Q[i,j] * (i == j ? x[i] : y[i,j]) for i=1:11, j=1:11) + Ξ²
)

@constraint(ising_ilp_model, c1[i=1:11,j=1:11;i!=j], y[i,j] >= x[i] + x[j] - 1)
@constraint(ising_ilp_model, c2[i=1:11,j=1:11;i!=j], y[i,j] <= x[i])
@constraint(ising_ilp_model, c3[i=1:11,j=1:11;i!=j], y[i,j] <= x[j])

println(ising_ilp_model)

Min -46 x[1] + 96 y[1,4] + 96 y[1,5] + 96 y[1,6] + 96 y[1,8] + 96 y[1,9] + 96 y[1,10] + 96 y[1,11] - 44 x[2] + 96 y[2,4] + 96 y[2,6] + 96 y[2,7] + 96 y[2,9] + 96 y[2,10] + 96 y[2,11] - 44 x[3] + 96 y[3,5] + 96 y[3,7] + 96 y[3,8] + 96 y[3,9] + 96 y[3,10] + 96 y[3,11] - 92 x[4] + 96 y[4,5] + 192 y[4,6] + 96 y[4,7] + 96 y[4,8] + 192 y[4,9] + 192 y[4,10] + 192 y[4,11] - 92 x[5] + 96 y[5,6] + 96 y[5,7] + 192 y[5,8] + 192 y[5,9] + 192 y[5,10] + 192 y[5,11] - 92 x[6] + 96 y[6,7] + 96 y[6,8] + 192 y[6,9] + 192 y[6,10] + 192 y[6,11] - 91 x[7] + 96 y[7,8] + 192 y[7,9] + 192 y[7,10] + 192 y[7,11] - 92 x[8] + 192 y[8,9] + 192 y[8,10] + 192 y[8,11] - 139 x[9] + 288 y[9,10] + 288 y[9,11] - 138 x[10] + 288 y[10,11] - 139 x[11] + 144

Subject to
 c1[1,2] : -x[1] - x[2] + y[1,2] β‰₯ -1
 c1[1,3] : -x[1] - x[3] + y[1,3] β‰₯ -1
 c1[1,4] : -x[1] - x[4] + y[1,4] β‰₯ -1
 c1[1,5] : -x[1] - x[5] + y[1,5] β‰₯ -1
 c1[1,6] : -x[1] - x[6] + y[1,6] β‰₯ -1
 c1[1,7] : -x[1] - x[7] + y[1,7] β‰₯ -1
 c1[1,8] : -x[1] - x[8] + y[1,8] β‰₯ -1
 c1[1,9] : -x[1] - x[9] + y[1,9] β‰₯ -1
 c1[1,10] : -x[1] - x[10] + y[1,10] β‰₯ -1
 c1[1,11] : -x[1] - x[11] + y[1,11] β‰₯ -1
 c1[2,1] : -x[1] - x[2] + y[2,1] β‰₯ -1
 c1[2,3] : -x[2] - x[3] + y[2,3] β‰₯ -1
 c1[2,4] : -x[2] - x[4] + y[2,4] β‰₯ -1
 c1[2,5] : -x[2] - x[5] + y[2,5] β‰₯ -1
 c1[2,6] : -x[2] - x[6] + y[2,6] β‰₯ -1
 c1[2,7] : -x[2] - x[7] + y[2,7] β‰₯ -1
 c1[2,8] : -x[2] - x[8] + y[2,8] β‰₯ -1
 c1[2,9] : -x[2] - x[9] + y[2,9] β‰₯ -1
 c1[2,10] : -x[2] - x[10] + y[2,10] β‰₯ -1
 c1[2,11] : -x[2] - x[11] + y[2,11] β‰₯ -1
 c1[3,1] : -x[1] - x[3] + y[3,1] β‰₯ -1
 c1[3,2] : -x[2] - x[3] + y[3,2] β‰₯ -1
 c1[3,4] : -x[3] - x[4] + y[3,4] β‰₯ -1
 c1[3,5] : -x[3] - x[5] + y[3,5] β‰₯ -1
 c1[3,6] : -x[3] - x[6] + y[3,6] β‰₯ -1
 c1[3,7] : -x[3] - x[7] + y[3,7] β‰₯ -1
 c1[3,8] : -x[3] - x[8] + y[3,8] β‰₯ -1
 c1[3,9] : -x[3] - x[9] + y[3,9] β‰₯ -1
 c1[3,10] : -x[3] - x[10] + y[3,10] β‰₯ -1
 c1[3,11] : -x[3] - x[11] + y[3,11] β‰₯ -1
 c1[4,1] : -x[1] - x[4] + y[4,1] β‰₯ -1
 c1[4,2] : -x[2] - x[4] + y[4,2] β‰₯ -1
 c1[4,3] : -x[3] - x[4] + y[4,3] β‰₯ -1
 c1[4,5] : -x[4] - x[5] + y[4,5] β‰₯ -1
 c1[4,6] : -x[4] - x[6] + y[4,6] β‰₯ -1
 c1[4,7] : -x[4] - x[7] + y[4,7] β‰₯ -1
 c1[4,8] : -x[4] - x[8] + y[4,8] β‰₯ -1
 c1[4,9] : -x[4] - x[9] + y[4,9] β‰₯ -1
 c1[4,10] : -x[4] - x[10] + y[4,10] β‰₯ -1
 c1[4,11] : -x[4] - x[11] + y[4,11] β‰₯ -1
 c1[5,1] : -x[1] - x[5] + y[5,1] β‰₯ -1
 c1[5,2] : -x[2] - x[5] + y[5,2] β‰₯ -1
 c1[5,3] : -x[3] - x[5] + y[5,3] β‰₯ -1
 c1[5,4] : -x[4] - x[5] + y[5,4] β‰₯ -1
 c1[5,6] : -x[5] - x[6] + y[5,6] β‰₯ -1
 c1[5,7] : -x[5] - x[7] + y[5,7] β‰₯ -1
 c1[5,8] : -x[5] - x[8] + y[5,8] β‰₯ -1
 c1[5,9] : -x[5] - x[9] + y[5,9] β‰₯ -1
 c1[5,10] : -x[5] - x[10] + y[5,10] β‰₯ -1
 c1[5,11] : -x[5] - x[11] + y[5,11] β‰₯ -1
[[...362 constraints skipped...]]
 y[6,7] binary
 y[7,7] binary
 y[8,7] binary
 y[9,7] binary
 y[10,7] binary
 y[11,7] binary
 y[1,8] binary
 y[2,8] binary
 y[3,8] binary
 y[4,8] binary
 y[5,8] binary
 y[6,8] binary
 y[7,8] binary
 y[8,8] binary
 y[9,8] binary
 y[10,8] binary
 y[11,8] binary
 y[1,9] binary
 y[2,9] binary
 y[3,9] binary
 y[4,9] binary
 y[5,9] binary
 y[6,9] binary
 y[7,9] binary
 y[8,9] binary
 y[9,9] binary
 y[10,9] binary
 y[11,9] binary
 y[1,10] binary
 y[2,10] binary
 y[3,10] binary
 y[4,10] binary
 y[5,10] binary
 y[6,10] binary
 y[7,10] binary
 y[8,10] binary
 y[9,10] binary
 y[10,10] binary
 y[11,10] binary
 y[1,11] binary
 y[2,11] binary
 y[3,11] binary
 y[4,11] binary
 y[5,11] binary
 y[6,11] binary
 y[7,11] binary
 y[8,11] binary
 y[9,11] binary
 y[10,11] binary
 y[11,11] binary


set_optimizer(ising_ilp_model, HiGHS.Optimizer)

optimize!(ising_ilp_model)

ising_ilp_x = round.(Int, value.(x))
ising_ilp_objective = objective_value(ising_ilp_model)

println(solution_summary(ising_ilp_model))
println("* objective value = $(round(ising_ilp_objective; digits=1))")
println("* x = $ising_ilp_x")
Running HiGHS 1.13.1 (git hash: 1d267d97c): Copyright (c) 2026 under Apache 2.0 license terms
Using BLAS: blastrampoline 
MIP has 330 rows; 132 cols; 770 nonzeros; 132 integer variables (132 binary)
Coefficient ranges:
  Matrix  [1e+00, 1e+00]
  Cost    [4e+01, 3e+02]
  Bound   [1e+00, 1e+00]
  RHS     [1e+00, 1e+00]
Presolving model
330 rows, 121 cols, 770 nonzeros  0s
165 rows, 66 cols, 385 nonzeros  0s
165 rows, 66 cols, 385 nonzeros  0s
Presolve reductions: rows 165(-165); columns 66(-66); nonzeros 385(-385) 
Objective function is integral with scale 1

Solving MIP model with:
   165 rows
   66 cols (66 binary, 0 integer, 0 implied int., 0 continuous, 0 domain fixed)
   385 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
     I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
     S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
     Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

 z       0       0         0   0.00%   -inf            144                Large        0      0      0         0     0.0s
         0       0         0   0.00%   -360.5          144              350.35%        0      0      5        13     0.0s
 L       0       0         0   0.00%   -140.9411765    5               2918.82%      228     18     10        73     0.1s

22.7% inactive integer columns, restarting
Model after restart has 135 rows, 51 cols (51 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 310 nonzeros

         0       0         0   0.00%   -124.8571429    5               2597.14%        6      0      0       118     0.1s
         0       0         0   0.00%   -124.8571429    5               2597.14%        6      6      2       131     0.1s

35.3% inactive integer columns, restarting
Model after restart has 99 rows, 33 cols (33 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 220 nonzeros

         0       0         0   0.00%   -86.5           5               1830.00%        8      0      0       157     0.1s
         0       0         0   0.00%   -86.5           5               1830.00%        8      8      5       170     0.1s

27.3% inactive integer columns, restarting
Model after restart has 60 rows, 24 cols (24 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 136 nonzeros

         0       0         0   0.00%   -77.16666667    5               1643.33%        5      0      0       183     0.1s
         0       0         0   0.00%   -77.16666667    5               1643.33%        5      5      3       199     0.1s
         1       0         1 100.00%   5               5                  0.00%       75      7      3       214     0.1s

Solving report
  Status            Optimal
  Primal bound      5
  Dual bound        5
  Gap               0% (tolerance: 0.01%)
  P-D integral      1.30449572572
  Solution status   feasible
                    5 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.11
  Max sub-MIP depth 2
  Nodes             1
  Repair LPs        0
  LP iterations     214
                    0 (strong br.)
                    125 (separation)
                    34 (heuristics)
solution_summary(; result = 1, verbose = false)

β”œ solver_name          : HiGHS
β”œ Termination
β”‚ β”œ termination_status : OPTIMAL
β”‚ β”œ result_count       : 1
β”‚ β”œ raw_status         : kHighsModelStatusOptimal
β”‚ β”” objective_bound    : 5.00000e+00
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”œ objective_value      : 5.00000e+00
β”‚ β”œ dual_objective_value : NaN
β”‚ β”” relative_gap         : 0.00000e+00
β”” Work counters
  β”œ solve_time (sec)   : 1.11032e-01
  β”œ simplex_iterations : 214
  β”œ barrier_iterations : -1
  β”” node_count         : 1
* objective value = 5.0
* x = [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]

Both the exact sampler and the HiGHS reformulation attain the minimum Ising objective value of 5. In spin variables, the two minimizing states have either s9=1s_9 = 1 or s11=1s_{11} = 1, with all other spins equal to -1. After converting the Ising model to the binary QUBO form for the integer-programming solve, those same minima appear as x9=1x_9 = 1 or x11=1x_{11} = 1, with all other binary variables equal to 0.

Let’s go back to the slidesΒΆ

Let’s now solve the QUBO problem using Simulated Annealing

set_optimizer(qubo_model, DWave.Neal.Optimizer)

set_optimizer_attribute(qubo_model, "num_reads", 1_000)

optimize!(qubo_model)

x = qubo_model[:x]

qubo_x = round.(Int, value.(x))
qubo_neal_objective = objective_value(qubo_model)

println(solution_summary(qubo_model))
println("* best sampled objective value = $(round(qubo_neal_objective; digits=1))")
println("* x = $qubo_x")

solution_summary(; result = 1, verbose = false)
β”œ solver_name          : D-Wave Neal Simulated Annealing Sampler
β”œ Termination
β”‚ β”œ termination_status : LOCALLY_SOLVED
β”‚ β”œ result_count       : 15
β”‚ β”” raw_status         : 
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”œ objective_value      : 5.00000e+00
β”‚ β”” dual_objective_value : 1.44000e+02
β”” Work counters
  β”” solve_time (sec)   : 4.35269e-01
* best sampled objective value = 5.0
* x = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]

plot(QUBOTools.EnergyFrequencyPlot(unsafe_backend(qubo_model)))
Loading...

Notice that this is the same example we have been solving earlier (via Integer Programming in the Quiz 1 and Ising model above).

Let’s go back to the slidesΒΆ

Let’s solve the graph coloring problem in the slides using QUBO.

Vertex kk-coloring of graphsΒΆ

Given a graph G(V,E)G(V, E), where VV is the set of vertices and EE is the set of edges of GG, and a positive integer kk, we ask if it is possible to assign a color to every vertex from VV, such that adjacent vertices have different colors assigned.

G(V,E)G(V, E) has 12 vertices and 23 edges. We ask if the graph is 3–colorable. Let’s first encode VV and EE using Julia’s built–in data structures:

Note: This second tutorial is heavily inspired by D-Wave’s Map coloring of Canada found here.

V = 1:12
E = [
    (1,2), (1,4), (1,6), (1,12),
    (2,3), (2,5), (2,7),
    (3,8), (3,10),
    (4,9), (4,11),
    (5,6), (5,9), (5,12),
    (6,7), (6,10),
    (7,8), (7,11),
    (8,9), (8,12),
    (9,10),
    (10,11),
    (11,12),
]

G = SimpleGraph(Edge.(E))
{12, 23} undirected simple Int64 graph
graph_layout = Vector{Point}(undef, 12)

graph_layout[1] = Point(-1.5,-1.5)
graph_layout[2] = Point(1.5,-1.5)
graph_layout[3] = Point(1.5,1.5)
graph_layout[4] = Point(-1.5,1.5)

for i in 5:12
    graph_layout[i] = Point(cos((2i+1) * Ο€/8), -sin((2i+1) * Ο€/8))
end

@drawsvg(
    begin
        sethue("black")
        background("white")
        drawgraph(
            G;
            layout=100 * graph_layout,
            vertexlabels = V,
        )
    end,
)
Loading...
# Valid configurations for the constraint that each node select a single color, in this case we want to use 3 colors
color_model = Model()

@variable(color_model, c[1:12,1:3], Bin)

# Each node must be colored with exactly one color
@constraint(color_model, unique[i=1:12], sum(c[i,:]) == 1)

# Add constraint that each pair of nodes with a shared edge not both select one color
@constraint(color_model, neigh[(i,j) ∈ E, k=1:3], c[i, k] * c[j,k] == 0)
2-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarQuadraticFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape},2,...} with index sets: Dimension 1, [(1, 2), (1, 4), (1, 6), (1, 12), (2, 3), (2, 5), (2, 7), (3, 8), (3, 10), (4, 9) … (5, 12), (6, 7), (6, 10), (7, 8), (7, 11), (8, 9), (8, 12), (9, 10), (10, 11), (11, 12)] Dimension 2, Base.OneTo(3) And data, a 23Γ—3 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarQuadraticFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}: neigh[(1, 2),1] : c[1,1]*c[2,1] = 0 … neigh[(1, 2),3] : c[1,3]*c[2,3] = 0 neigh[(1, 4),1] : c[1,1]*c[4,1] = 0 neigh[(1, 4),3] : c[1,3]*c[4,3] = 0 neigh[(1, 6),1] : c[1,1]*c[6,1] = 0 neigh[(1, 6),3] : c[1,3]*c[6,3] = 0 neigh[(1, 12),1] : c[1,1]*c[12,1] = 0 neigh[(1, 12),3] : c[1,3]*c[12,3] = 0 neigh[(2, 3),1] : c[2,1]*c[3,1] = 0 neigh[(2, 3),3] : c[2,3]*c[3,3] = 0 neigh[(2, 5),1] : c[2,1]*c[5,1] = 0 … neigh[(2, 5),3] : c[2,3]*c[5,3] = 0 neigh[(2, 7),1] : c[2,1]*c[7,1] = 0 neigh[(2, 7),3] : c[2,3]*c[7,3] = 0 neigh[(3, 8),1] : c[3,1]*c[8,1] = 0 neigh[(3, 8),3] : c[3,3]*c[8,3] = 0 neigh[(3, 10),1] : c[3,1]*c[10,1] = 0 neigh[(3, 10),3] : c[3,3]*c[10,3] = 0 neigh[(4, 9),1] : c[4,1]*c[9,1] = 0 neigh[(4, 9),3] : c[4,3]*c[9,3] = 0 neigh[(4, 11),1] : c[4,1]*c[11,1] = 0 … neigh[(4, 11),3] : c[4,3]*c[11,3] = 0 neigh[(5, 6),1] : c[5,1]*c[6,1] = 0 neigh[(5, 6),3] : c[5,3]*c[6,3] = 0 neigh[(5, 9),1] : c[5,1]*c[9,1] = 0 neigh[(5, 9),3] : c[5,3]*c[9,3] = 0 neigh[(5, 12),1] : c[5,1]*c[12,1] = 0 neigh[(5, 12),3] : c[5,3]*c[12,3] = 0 neigh[(6, 7),1] : c[6,1]*c[7,1] = 0 neigh[(6, 7),3] : c[6,3]*c[7,3] = 0 neigh[(6, 10),1] : c[6,1]*c[10,1] = 0 … neigh[(6, 10),3] : c[6,3]*c[10,3] = 0 neigh[(7, 8),1] : c[7,1]*c[8,1] = 0 neigh[(7, 8),3] : c[7,3]*c[8,3] = 0 neigh[(7, 11),1] : c[7,1]*c[11,1] = 0 neigh[(7, 11),3] : c[7,3]*c[11,3] = 0 neigh[(8, 9),1] : c[8,1]*c[9,1] = 0 neigh[(8, 9),3] : c[8,3]*c[9,3] = 0 neigh[(8, 12),1] : c[8,1]*c[12,1] = 0 neigh[(8, 12),3] : c[8,3]*c[12,3] = 0 neigh[(9, 10),1] : c[9,1]*c[10,1] = 0 … neigh[(9, 10),3] : c[9,3]*c[10,3] = 0 neigh[(10, 11),1] : c[10,1]*c[11,1] = 0 neigh[(10, 11),3] : c[10,3]*c[11,3] = 0 neigh[(11, 12),1] : c[11,1]*c[12,1] = 0 neigh[(11, 12),3] : c[11,3]*c[12,3] = 0
set_optimizer(color_model, () -> ToQUBO.Optimizer(DWave.Neal.Optimizer))

optimize!(color_model)

color_ρ = get_optimizer_attribute.(neigh, ToQUBO.Attributes.ConstraintEncodingPenalty())

color_c = round.(Int, value.(c))

println(solution_summary(color_model))
println("* c = $color_c")
println("* ρ = $color_ρ")

solution_summary(; result = 1, verbose = false)
β”œ solver_name          : Virtual QUBO Model
β”œ Termination
β”‚ β”œ termination_status : LOCALLY_SOLVED
β”‚ β”œ result_count       : 170
β”‚ β”” raw_status         : 
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : FEASIBLE_POINT
β”‚ β”œ dual_status          : NO_SOLUTION
β”‚ β”” objective_value      : 0.00000e+00
β”” Work counters
  β”” solve_time (sec)   : 5.54342e-01
* c = [1 0 0; 0 0 1; 0 1 0; 0 1 0; 1 0 0; 0 0 1; 0 1 0; 1 0 0; 0 0 1; 1 0 0; 0 0 1; 0 1 0]
* ρ = 2-dimensional DenseAxisArray{Float64,2,...} with index sets:
    Dimension 1, [(1, 2), (1, 4), (1, 6), (1, 12), (2, 3), (2, 5), (2, 7), (3, 8), (3, 10), (4, 9)  …  (5, 12), (6, 7), (6, 10), (7, 8), (7, 11), (8, 9), (8, 12), (9, 10), (10, 11), (11, 12)]
    Dimension 2, Base.OneTo(3)
And data, a 23Γ—3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0


plot(QUBOTools.EnergyFrequencyPlot(unsafe_backend(color_model).optimizer))
Loading...
@drawsvg(
    begin
        sethue("black")
        background("white")
        drawgraph(
            G;
            layout=100 * graph_layout,
            vertexlabels = V,
            vertexfillcolors = (v) -> begin
                if color_c[v,1] > 0
                    return colorant"red"
                elseif color_c[v,2] > 0
                    return colorant"blue"
                elseif color_c[v,3] > 0
                    return colorant"green"
                else
                    return colorant"white"
                end
            end
        )
    end,
)
Loading...
QUBO
QUBO and Ising Models (Python)
GAMA
Graver Augmentation Multiseed Algorithm (Python)