Experimental MILP Solvers — SCIP, HiGHS, GLPK, CBC

PyQBPP can solve QUBO expressions with several third-party exact MILP solvers. They share a single interface, so user code can switch between them — and with pyqbpp.ABS3Solver — by changing only the class name.

These solvers minimize a linear objective, so the quadratic QUBO must be linearized before it is handed over (see below). Solvers that accept the quadratic objective directly do not appear here; they are documented under QUBO/HUBO Solvers. The constraint-programming engine OR-Tools CP-SAT is documented under CP Solvers.

Experimental. These integrations are provided for experimentation and benchmarking. Their API may change without notice, and each solver’s Python binding must be installed separately (see Setup). They solve QUBO (degree ≤ 2) only — reduce a HUBO to QUBO first, or use pyqbpp.ABS3Solver / pyqbpp.EasySolver, which support arbitrary degree.

Internally each binary product is replaced by an auxiliary variable with linear linking constraints (Fortet linearization), so the QUBO is handed to the solver as a pure MILP. The returned solution’s energy is always recomputed exactly from the original QUBO.

Solver	Class	Python binding	License
SCIP	`pyqbpp.ScipSolver`	PySCIPOpt	Apache-2.0
HiGHS	`pyqbpp.HighsSolver`	highspy	MIT
GLPK	`pyqbpp.GlpkSolver`	swiglpk	GPL
CBC	`pyqbpp.CbcSolver`	python-mip	EPL

For the exact solvers that accept the quadratic objective directly (DimodExactSolver and QiskitOptimizationSolver), see QUBO/HUBO Solvers.

Usage

All four solvers share the same interface. The following program solves a number partitioning problem with SCIP — replace qbpp.ScipSolver with qbpp.HighsSolver, qbpp.GlpkSolver, or qbpp.CbcSolver to use a different one:

import pyqbpp as qbpp

w = qbpp.array([64, 27, 47, 74, 12, 83, 63, 40])
x = qbpp.var("x", shape=len(w))
p = 0
q = 0
for i in range(len(w)):
    p += w[i] * x[i]
    q += w[i] * ~x[i]
f = qbpp.sqr(p - q)
f.simplify_as_binary()

solver = qbpp.ScipSolver(f)            # swap the class name to switch solver
sol = solver.search(time_limit=10.0)

print(f"energy = {sol.energy}")
print(f"bound  = {sol.info.get('bound')}")
print(f"status = {sol.info.get('status')}")
print("P:", [w[i] for i in range(len(w)) if sol(x[i]) == 1])
print("Q:", [w[i] for i in range(len(w)) if sol(x[i]) == 0])

When the energy equals the bound, the solution is proven optimal. A solver object is created from an expression; on construction the QUBO is linearized into the solver’s internal MILP model. A higher-degree (HUBO) expression raises an exception.

ILP mode

A model whose objective is linear and whose constraints are written with qbpp.cons() (it may contain native integer variables) can be handed over as an integer linear program, bypassing the QUBO linearization, by passing ilp=True (qbpp.ScipSolver(f, ilp=True) — the same for Highs / Cbc / Glpk). Integer variables become integer columns and qbpp.cons() constraints become linear rows (weights are ignored; the constraints are hard). The constructor raises RuntimeError when the objective has a term of degree 2 or higher, when a constraint body is nonlinear (degree 2 or higher), or when an allowed-value set equal=[...] is used. Without ilp=True (the plain mode), qbpp.cons() and native integer variables are not accepted (QUBO linearization only).

Parameters

Parameters are passed to search() as keyword arguments (a dict is also accepted). The keys accepted by every wrapper are (a key marked with solvers in parentheses is ignored by the others):

Key	Value	Description
`time_limit`	seconds	Stop when the time limit is reached
`target_energy`	energy	Stop when an energy ≤ this value is found
`callback_timer_interval`	seconds	Initial interval for `Timer` events
`enable_default_callback`	`1`	Print energy and TTS of each new incumbent (SCIP/HiGHS)
`thread_count`	threads	Worker threads (SCIP/HiGHS/CBC)
`topk_sols`	K	Return up to K solutions (best-effort, SCIP/HiGHS)
`gap_limit`	gap	Relative MIP gap stop (SCIP/HiGHS)
`output_flag`	`1`	Show the solver’s own log (SCIP/HiGHS)

Solver-specific extras:

SCIP — any unrecognized key is forwarded to SCIP verbatim (e.g. solver.search({"limits/gap": 0.0})). The formulation (linearize / quadratic) is a constructor argument, not a search() key — see below.
HiGHS — any unrecognized key is forwarded to HiGHS setOptionValue (e.g. presolve="on").

SCIP: linearize vs quadratic formulation

ScipSolver can feed the QUBO to SCIP in two ways:

"linearize" (default) — Fortet linearization to a pure MILP (the same transform used by the other solvers), giving SCIP a tight LP relaxation.
"quadratic" — SCIP’s objective is strictly linear, so the QUBO is modelled with one extra objective variable and a single quadratic (nonlinear) constraint that SCIP’s nonlinear handler reformulates itself, adding no Fortet auxiliary variables. Its per-term relaxation is weaker, so it is usually slower on dense penalty QUBOs — provided mainly for comparison.

Both reach the same optimum. The formulation is fixed at construction — to use a different one, create a new solver object (it is not a search() keyword):

solver = qbpp.ScipSolver(f, formulation="quadratic")
sol = solver.search()

This option is specific to SCIP; HiGHS, GLPK, and CBC always use the linearized MILP.

Solver Info

sol.info carries solver-produced strings: status, bound, mip_gap (SCIP/HiGHS), node_count, solution_count, <solver>_version (scip_version / highs_version / glpk_version), and run_time. solution_count is the number of solutions the solver holds for SCIP and CBC, and 1 if a solution was found (else 0) for HiGHS and GLPK.

Custom Callback

The callback API is identical to qbpp.ABS3Solver, and there are two ways to use it:

Register a function — solver.set_callback(fn). fn is called as fn(solver), receiving the solver itself (no inheritance)
Subclass — derive from the solver and override callback()

If both are used, the override wins. Call super().callback() from the override to run both.

Either way, event(), best_sol() (a pyqbpp.Sol), bound(), timer(seconds), and terminate() are available. The event constants are class attributes: EVENT_START, EVENT_BEST_UPDATED, EVENT_TIMER.

Method 1: registering a function

import pyqbpp as qbpp

def on_event(s):
    if s.event() == s.EVENT_BEST_UPDATED:
        best = s.best_sol()
        print(f"New best: energy={best.energy} TTS={best.tts}s bound={s.bound()}")
        if best.energy == 0:
            s.terminate()             # stop as soon as the optimum is found

x = qbpp.var("x", shape=8)
f = qbpp.sqr(qbpp.sum(x) - 4)
f.simplify_as_binary()
solver = qbpp.HighsSolver(f)
solver.set_callback(on_event)
sol = solver.search(time_limit=5)
print(f"energy={sol.energy}")

Method 2: subclassing

import pyqbpp as qbpp

class MySolver(qbpp.HighsSolver):
    def callback(self):
        if self.event() == self.EVENT_BEST_UPDATED:
            s = self.best_sol()
            print(f"New best: energy={s.energy} TTS={s.tts}s bound={self.bound()}")
            if s.energy == 0:
                self.terminate()          # stop as soon as optimal is found

x = qbpp.var("x", shape=8)
f = qbpp.sqr(qbpp.sum(x) - 4)
f.simplify_as_binary()
sol = MySolver(f).search(time_limit=5)
print(f"energy={sol.energy}")

Live-callback support differs by binding. ScipSolver and HighsSolver deliver BestUpdated/Timer during the solve. GlpkSolver and CbcSolver deliver only EVENT_START — swiglpk cannot install a Python branch-and-cut callback and python-mip’s incumbent callback is not invoked on common builds. The C++ qbpp::GlpkSolver / qbpp::CbcSolver do deliver live events. hint(sol) warm-starts SCIP/HiGHS; it is a no-op for GLPK/CBC. Use time_limit to bound a run.

Setup

Install only the bindings you use:

Solver	Install
SCIP	`conda install -c conda-forge pyscipopt`
HiGHS	`pip install highspy` (or `conda install -c conda-forge highspy`)
GLPK	`conda install -c conda-forge swiglpk`
CBC	`pip install mip`

The bindings are imported lazily, so import pyqbpp works even when a solver is not installed; the error is raised only when you construct that solver.