実験的 MILP ソルバー — SCIP, HiGHS, GLPK, CBC¶

Hi-QUBO は複数のサードパーティ製厳密 MILP ソルバーで QUBO 式を解くことができます.これらは共通インタフェースを持つヘッダオンリーのソルバーとしてラップされており,クラス名を変えるだけで互いに,また qbpp::ABS3Solver とも切り替えられます.

これらは線形の目的関数を最小化するため,二次の QUBO は渡す前に 線形化する必要があります（後述）.これがこのページの判定基準です.

制約プログラミングエンジンの OR-Tools CP-SAT は CP ソルバーを参照してください.

実験的機能. これらは実験・ベンチマーク用途で提供されます.API は予告なく変更される可能性があり,各ソルバーは別途インストールが必要です（セットアップ参照）. 対応は QUBO（次数 ≤ 2）のみです.HUBO は事前に QUBO へ削減するか,任意次数に対応する qbpp::ABS3Solver / qbpp::EasySolver を使ってください.

内部では各二次項 x·y を補助変数＋線形リンク制約に置き換え（Fortet 線形化）, QUBO を純粋な MILP としてソルバーに渡します.返される解のエネルギーは,ソルバーの浮動小数点目的値とは独立に,元の QUBO から常に厳密に再計算されます.

ソルバー	クラス	ライセンス	備考
SCIP	`qbpp::ScipSolver`	Apache-2.0 (OSS)	linearize / quadratic 定式化
HiGHS	`qbpp::HighsSolver`	MIT (OSS)	高速な OSS MILP
GLPK	`qbpp::GlpkSolver`	GPL (OSS)	軽量
CBC	`qbpp::CbcSolver`	EPL (OSS)	COIN-OR branch & cut

使い方

4 つのソルバーはすべて同じインタフェースです.次のプログラムは SCIP で数分割問題を解きます.qbpp::ScipSolver を qbpp::HighsSolver / qbpp::GlpkSolver / qbpp::CbcSolver に置き換える（と対応ヘッダを include する）だけで別のソルバーを使えます:

milp-solvers-program1.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/scip.hpp>     // または highs.hpp / glpk.hpp / cbc.hpp

int main() {
  std::vector<int> w = {64, 27, 47, 74, 12, 83, 63, 40};
  auto x = qbpp::var("x", w.size());
  auto p = qbpp::toExpr(0), q = qbpp::toExpr(0);
  for (size_t i = 0; i < w.size(); ++i) { p += w[i] * x[i]; q += w[i] * ~x[i]; }
  auto f = qbpp::sqr(p - q);
  f.simplify_as_binary();

  auto solver = qbpp::ScipSolver(f);   // <- クラス名を差し替えればソルバー切替
  auto sol = solver.search({{"time_limit", 10.0}});

  std::cout << "energy = " << sol.energy() << std::endl;
  std::cout << "bound  = " << sol.info("bound") << std::endl;
  std::cout << "status = " << sol.info("status") << std::endl;
}

エネルギーが下界 sol.info("bound") と一致すれば,その解は最適であることが保証されます:

energy = 0
bound  = 0.000000
status = OPTIMAL

ソルバーオブジェクトは式から生成され,構築時に（simplify 済みの）QUBO がソルバー内部の MILP モデルへ線形化されます.式が高次（HUBO）項を含む場合,コンストラクタは例外を送出します.

パラメータ

パラメータは search() にキー／値のペアの初期化子リストで渡します.すべてのラッパーが解釈する共通キーは次の通りです:

キー	値	説明
`time_limit`	秒	制限時間に達したら停止
`target_energy`	エネルギー	この値以下の解を見つけたら停止
`callback_timer_interval`	秒	`Timer` イベントの初期間隔
`enable_default_callback`	`1`	新しい incumbent ごとにエネルギーと TTS を表示
`thread_count`	スレッド数	ワーカースレッド数（SCIP/HiGHS.GLPK/CBC は無視）
`topk_sols`	K	最大 K 個の解を返す（ベストエフォート,下記参照）
`gap_limit`	gap	相対 MIP gap による停止（SCIP/HiGHS）
`output_flag`	`1`	ソルバー自身のログを表示（SCIP/HiGHS）

ソルバー固有の追加:

SCIP — 未知のキーは SCIP にそのまま転送されます（例 "limits/gap", "lp/threads"）.formulation（linearize / quadratic）は コンストラクタのオプションであり,search() のキーではありません（下記参照）.
HiGHS — 未知のキーは HiGHS の setOptionValue に転送されます（例 "presolve", "mip_rel_gap"）.

topk_sols はベストエフォートです.これらのソルバーは内部ストレージに残る相異なる解（多くの場合 incumbent のみ）を返します.

SCIP: linearize と quadratic 定式化

qbpp::ScipSolver は QUBO を 2 通りの方法で SCIP に渡せます:

"linearize"（既定） — Fortet 線形化で純粋な MILP にする（他のソルバーと共通の変換）.SCIP に締まった LP 緩和を与えます.
"quadratic" — SCIP の目的関数は線形のみのため,目的変数 t を 1 つ追加し, 2 次（非線形）制約 t == const + Σ qᵢⱼ·xᵢ·xⱼ を 1 本張って t を最小化します. 二次項は SCIP の非線形制約ハンドラが内部で再定式化するため,Fortet 補助変数は追加されません.ただし項ごとの（McCormick）緩和は緩く,密なペナルティ QUBO では通常遅くなります.比較用に提供しています.

どちらの定式化でも同じ最適解に到達します.定式化は構築時に固定されます. 別の定式化を使うにはソルバーオブジェクトを作り直してください（search() のパラメータではありません）:

qbpp::ScipSolver solver(f, qbpp::ScipSolver::Formulation::Quadratic);
auto sol = solver.search();

このオプションは SCIP 固有です.HiGHS・GLPK・CBC は常に線形化 MILP を使います.

Solver Info

sol.info() はソルバーが生成した文字列を保持します:

キー	説明
`status`	`OPTIMAL`, `TIME_LIMIT`, `INFEASIBLE`, …（綴りはソルバー依存）
`bound`	最良の双対下界
`mip_gap`	最終的な相対 MIP gap（SCIP/HiGHS）
`node_count`	分枝限定ノード数
`solution_count`	解が得られたら `1`,なければ `0`
`<solver>_version`	`scip_version` / `highs_version` / `glpk_version`
`run_time`	実時間の求解時間（秒）

カスタムコールバック

コールバック API は qbpp::ABS3Solver と同一で、次の 2 つの方法があります。

関数を登録する — solver.set_callback(fn)。fn はソルバー自身を引数として呼ばれます（継承不要）
サブクラス化する — ソルバーを継承し callback() 仮想メソッドを override します

両方を指定した場合はオーバーライドが優先されます。両方実行したい場合はオーバーライドの中で基底クラスの callback() を呼んでください。

イベントは以下の 3 種です:

イベント	説明
`CallbackEvent::Start`	`search()` の開始時に 1 度呼ばれる
`CallbackEvent::BestUpdated`	新しい incumbent が見つかるたびに呼ばれる
`CallbackEvent::Timer`	`timer(seconds)` で設定した間隔で定期的に呼ばれる

コールバック内では event(),best_sol()（現在の最良 qbpp::Sol）, bound()（現在の双対下界）,timer(seconds)（タイマー設定/無効化）, terminate()（次の安全点で探索を停止）が使えます.hint(sol) は SCIP と HiGHS をウォームスタートします（GLPK・CBC では何もしません）.

方法1: 関数を登録する

milp-solvers-program2.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/highs.hpp>

int main() {
  auto x = qbpp::var("x", 8);
  auto f = qbpp::sqr(qbpp::sum(x) - 4);
  f.simplify_as_binary();

  qbpp::HighsSolver solver(f);
  solver.set_callback([](const qbpp::HighsSolver& s) {
    if (s.event() == qbpp::CallbackEvent::BestUpdated) {
      std::cout << "New best: energy=" << s.best_sol().energy()
                << "  bound=" << s.bound() << std::endl;
      if (s.best_sol().energy() == 0) s.terminate();
    }
  });

  auto sol = solver.search({{"time_limit", 5}});
  std::cout << "energy=" << sol.energy() << std::endl;
}

方法2: サブクラス化する

milp-solvers-program3.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/highs.hpp>

class MySolver : public qbpp::HighsSolver {
 public:
  using HighsSolver::HighsSolver;
  void callback() const override {
    if (event() == qbpp::CallbackEvent::BestUpdated) {
      std::cout << "New best: energy=" << best_sol().energy()
                << " TTS=" << best_sol().tts() << "s  bound=" << bound() << std::endl;
      if (best_sol().energy() == 0) terminate();  // 最適が見つかり次第停止
    }
  }
};

int main() {
  auto x = qbpp::var("x", 8);
  auto f = qbpp::sqr(qbpp::sum(x) - 4);
  f.simplify_as_binary();
  auto sol = MySolver(f).search({{"time_limit", 5}});
  std::cout << "energy=" << sol.energy() << std::endl;
}

セットアップ

各ソルバーは別途インストールが必要です.ヘッダは独立しています（<qbpp/scip.hpp>, <qbpp/highs.hpp>, <qbpp/glpk.hpp>, <qbpp/cbc.hpp>）. 使うものだけ include してください.qbpp 本体は dlopen で読み込まれるため -lqbpp は不要です.4 つすべてを手軽に入れるには conda-forge が便利です:

conda install -c conda-forge scip highs glpk coincbc

ソルバー別のビルドフラグ（他の qbpp プログラム同様 -ldl -pthread を付与）:

ソルバー	リンクフラグ	ヘッダパス（conda）
SCIP	`-lscip`	（システム include か `-I$PREFIX/include`）
HiGHS	`-lhighs`	`-isystem $PREFIX/include/highs`
GLPK	`-lglpk`	（システム include か `-I$PREFIX/include`）
CBC	`-lCbc -lCbcSolver -lCgl -lOsiClp -lClp -lOsi -lCoinUtils`	`-isystem $PREFIX/include/coin`

例（conda を $CONDA_PREFIX に導入した場合）:

g++ -std=c++17 -Wfloat-conversion your_program.cpp -o your_program \
    -isystem $CONDA_PREFIX/include/highs \
    -L$CONDA_PREFIX/lib -Wl,-rpath,$CONDA_PREFIX/lib -lhighs -ldl -pthread

SCIP は apt/deb（SCIPOptSuite-*.deb）で導入するとヘッダと libscip.so が既定のパスに置かれるため,-lscip だけで済みます.

PyQBPP は複数のサードパーティ製厳密 MILP ソルバーで QUBO 式を解くことができます.これらは共通インタフェースを持ち,クラス名を変えるだけで互いに, また pyqbpp.ABS3Solver とも切り替えられます.

これらは線形の目的関数を最小化するため,二次の QUBO は渡す前に 線形化する必要があります（後述）.二次目的関数を直接受け取れるソルバーはここには含まれず, QUBO/HUBO ソルバーにまとめています.制約プログラミングエンジンの OR-Tools CP-SAT は CP ソルバーを参照してください.

実験的機能. これらは実験・ベンチマーク用途で提供されます.API は予告なく変更される可能性があり,各ソルバーの Python バインディングは別途インストールが必要です（セットアップ参照）.対応は QUBO（次数 ≤ 2）のみです. HUBO は事前に QUBO へ削減するか,任意次数に対応する pyqbpp.ABS3Solver / pyqbpp.EasySolver を使ってください.

内部では各二次項を補助変数＋線形リンク制約に置き換え（Fortet 線形化）,QUBO を純粋な MILP としてソルバーに渡します.返される解のエネルギーは元の QUBO から常に厳密に再計算されます.

ソルバー	クラス	Python バインディング	ライセンス
SCIP	`pyqbpp.ScipSolver`	PySCIPOpt	Apache-2.0
HiGHS	`pyqbpp.HighsSolver`	highspy	MIT
GLPK	`pyqbpp.GlpkSolver`	swiglpk	GPL
CBC	`pyqbpp.CbcSolver`	python-mip	EPL

二次目的関数を直接受け取れる商用の厳密ソルバーは QUBO/HUBO ソルバーを参照してください.

使い方

4 つのソルバーはすべて同じインタフェースです.次のプログラムは SCIP で数分割問題を解きます.qbpp.ScipSolver を qbpp.HighsSolver / qbpp.GlpkSolver / qbpp.CbcSolver に置き換えるだけで別のソルバーを使えます:

milp-solvers-program1.py¶

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)            # クラス名を差し替えればソルバー切替
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])

エネルギーが下界と一致すれば最適が保証されます.ソルバーオブジェクトは式から生成され,構築時に QUBO がソルバー内部の MILP モデルへ線形化されます.高次（HUBO）の式は例外を送出します.

パラメータ

パラメータは search() にキーワード引数で渡します（dict も可）.すべてのラッパーが解釈する共通キーは次の通りです:

キー	値	説明
`time_limit`	秒	制限時間に達したら停止
`target_energy`	エネルギー	この値以下の解を見つけたら停止
`callback_timer_interval`	秒	`Timer` イベントの初期間隔
`enable_default_callback`	`1`	新しい incumbent ごとにエネルギーと TTS を表示
`thread_count`	スレッド数	ワーカースレッド数（SCIP/HiGHS/CBC）
`topk_sols`	K	最大 K 個の解を返す（ベストエフォート）
`gap_limit`	gap	相対 MIP gap による停止（SCIP/HiGHS）
`output_flag`	`1`	ソルバー自身のログを表示（SCIP/HiGHS）

ソルバー固有の追加:

SCIP — 未知のキーは SCIP にそのまま転送されます（例 solver.search({"limits/gap": 0.0})）.formulation（linearize / quadratic）は コンストラクタの引数であり,search() のキーではありません（下記参照）.
HiGHS — 未知のキーは HiGHS の setOptionValue に転送されます（例 presolve="on"）.

SCIP: linearize と quadratic 定式化

ScipSolver は QUBO を 2 通りの方法で SCIP に渡せます:

"linearize"（既定） — Fortet 線形化で純粋な MILP にする（他のソルバーと共通の変換）.SCIP に締まった LP 緩和を与えます.
"quadratic" — SCIP の目的関数は線形のみのため,目的変数を 1 つ追加し,SCIP の非線形制約ハンドラが内部で再定式化する 2 次（非線形）制約を 1 本張ります. Fortet 補助変数は追加されません.項ごとの緩和は緩いため,密なペナルティ QUBO では通常遅くなります.比較用に提供しています.

どちらも同じ最適解に到達します.定式化は構築時に固定されます.別の定式化を使うにはソルバーオブジェクトを作り直してください（search() のキーワードではありません）:

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

このオプションは SCIP 固有です.HiGHS・GLPK・CBC は常に線形化 MILP を使います.

Solver Info

sol.info はソルバーが生成した文字列を保持します: status, bound, mip_gap （SCIP/HiGHS）, node_count, solution_count, <solver>_version （scip_version / highs_version / glpk_version）, run_time.

カスタムコールバック

コールバック API は qbpp.ABS3Solver と同一で、次の 2 つの方法があります。

関数を登録する — solver.set_callback(fn)。fn はソルバー自身を引数として fn(solver) の形で呼ばれます（継承不要）
サブクラス化する — ソルバーを継承し callback() を override します

両方を指定した場合はオーバーライドが優先されます。両方実行したい場合はオーバーライドの中で super().callback() を呼んでください。

どちらの場合も event()、best_sol()（pyqbpp.Sol）、bound()、timer(seconds)、 terminate() が使えます。イベント定数はクラス属性 EVENT_START / EVENT_BEST_UPDATED / EVENT_TIMER です.

方法1: 関数を登録する

milp-solvers-program2.py¶

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()             # 最適が見つかり次第停止

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}")

方法2: サブクラス化する

milp-solvers-program3.py¶

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()          # 最適が見つかり次第停止

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}")

ライブコールバックの対応はバインディング依存です. ScipSolver と HighsSolver は求解中に BestUpdated/Timer を発火します.GlpkSolver と CbcSolver は EVENT_START のみです（swiglpk は Python コールバックを設定できず,python-mip の incumbent コールバックも一般的なビルドでは呼ばれないため）. C++ の qbpp::GlpkSolver / qbpp::CbcSolver はライブイベントを発火します. hint(sol) は SCIP/HiGHS をウォームスタートします（GLPK/CBC では何もしません）. 求解時間は time_limit で制限してください.

セットアップ

使うバインディングだけ入れてください:

ソルバー	インストール
SCIP	`conda install -c conda-forge pyscipopt`
HiGHS	`pip install highspy`（または `conda install -c conda-forge highspy`）
GLPK	`conda install -c conda-forge swiglpk`
CBC	`pip install mip`

バインディングは遅延 import されるため,未インストールでも import pyqbpp は成功します.エラーはそのソルバーを生成した時点で初めて送出されます.