ネイティブ制約¶
Hi-QUBO では,式の中の制約部分を qbpp::cons() で囲むと,その部分は
制約とみなされて特別に処理されます.Hi-QUBO にバンドルされている
ソルバーは,宣言された制約を満たすように効率よく探索を行います.
整数線形計画法を cons() で解く
範囲制約と整数線形計画法の求解では,次の整数線形計画問題を,
範囲制約 c1・c2 を重み付きのペナルティ式として目的関数に加える方法で
解きました:
同じ問題は,制約を qbpp::cons() で囲むと次のように書けます:
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>
int main() {
auto x = 0 <= qbpp::var_int("x") <= 10;
auto y = 0 <= qbpp::var_int("y") <= 10;
auto f = 5 * x + 4 * y;
auto c1 = 0 <= 2 * x + 3 * y <= 24;
auto c2 = 0 <= 7 * x + 5 * y <= 54;
auto g = -f + 100 * (qbpp::cons(c1) + qbpp::cons(c2));
g.simplify_as_binary();
auto solver = qbpp::EasySolver(g);
auto sol = solver.search({{"time_limit", 1.0}});
std::cout << "x = " << sol(x) << ", y = " << sol(y) << std::endl;
std::cout << "f = " << sol(f) << std::endl;
std::cout << "violated constraints = " << g.cons(sol) << std::endl;
}
変更点は,ペナルティ和 100 * (c1 + c2) を
100 * (qbpp::cons(c1) + qbpp::cons(c2)) に書き換えただけです.
これだけで c1・c2 は単なるペナルティ式ではなく制約として宣言され,
ソルバーは制約を満たす解を効率よく探索します.
g.cons(sol) は解 sol で違反している制約の本数を返します(0 なら全制約を充足).
プログラムの出力は以下の通りです:
x = 4, y = 5
f = 40
violated constraints = 0
ナップサック問題の例
もう1つの例として,次のプログラムは簡単なナップサック問題
(容量制約と等式制約)を cons() で解きます:
#include <qbpp/easy_solver.hpp>
#include <qbpp/qbpp.hpp>
int main() {
auto x = qbpp::var("x", 6);
int value[] = {3, 5, 2, 7, 4, 6};
int weight[] = {2, 4, 1, 5, 3, 4};
qbpp::Expr obj, load;
for (int i = 0; i < 6; ++i) {
obj += -value[i] * x(i);
load += weight[i] * x(i);
}
auto f = obj + 100 * qbpp::cons(load <= 8) // capacity, weight 100
+ 10 * qbpp::cons(x(0) + x(1) == 1); // equality, weight 10
f.simplify_as_binary();
qbpp::EasySolver solver(f);
auto sol = solver.search({{"time_limit", 1}});
std::cout << "objective = " << sol.energy() << std::endl;
}
従来のペナルティ式からの移行は,制約部分を qbpp::cons() で囲むだけです —
obj + 1000 * (rows + cols) を obj + 1000 * qbpp::cons(rows + cols) に
書き換えます.多くの問題で,同じ制約をペナルティ式のまま解くより
大幅に良い解が得られます.
制約の書き方
制約は右辺が整数の比較,または連鎖した両側範囲を qbpp::cons() で囲んで
書きます.重みは制約へのスカラー係数として書き,+ で目的関数や
他の制約と自由に足し合わせられます.
auto a = qbpp::var("a"), b = qbpp::var("b"), c = qbpp::var("c");
auto u = qbpp::var("u"), v = qbpp::var("v"), w = qbpp::var("w");
auto cons = 1000 * qbpp::cons(a + b + c == 2) // equality
+ 1000 * qbpp::cons(10 * u + 30 * v - 10 * w <= 35) // one-sided
+ 500 * qbpp::cons(0 <= a + b - c <= 1); // two-sided
配列の比較を囲むと要素ごとに 1 本の制約になるので,行列の one-hot 行は 1 文で書けます.
auto y = qbpp::var("y", 4, 4);
auto one_hot = 1000 * qbpp::cons(qbpp::vector_sum(y) == 1); // one per row
制約の重み付き和をまとめて qbpp::cons() で囲むこともできます.次の例は
重み 100 と 150 の 2 本の制約を作ります.
auto p = qbpp::var("p"), q = qbpp::var("q"), r = qbpp::var("r");
auto cons2 = qbpp::cons(100 * (p + q + r == 2) +
150 * (10 * p + 30 * q - 10 * r <= 35));
制約を蓄積した式に *= を使うと,蓄積済みの全制約の重みを一括で
スケールできます.
auto k0 = qbpp::var("k0"), k1 = qbpp::var("k1"), k2 = qbpp::var("k2");
auto cons4 = qbpp::cons(k0 + k1 + k2 == 2);
cons4 += qbpp::cons(k0 + k1 - k2 <= 1);
cons4 *= 1000; // scale ALL weights at once
式を出力すると目的関数の多項式が,f.cons() を出力すると宣言された
制約リストが表示されます(重みが 1 のときは係数プレフィックスを省略,
片側制約は片側表示).
auto m0 = qbpp::var("m0"), m1 = qbpp::var("m1"), m2 = qbpp::var("m2");
auto printed = 1000 * qbpp::cons(m0 + m1 + m2 == 2)
+ 500 * qbpp::cons(0 <= m0 + m1 - m2 <= 1);
std::cout << printed.cons() << std::endl;
出力は次のようになります.
1000 * (m0 +m1 +m2 == 2)
500 * (0 <= m0 +m1 -m2 <= 1)
離散許容値集合
式の値がとびとびの許容値のいずれかに一致することを要求する制約は,
qbpp::equal{...} で書けます.qbpp::cons(s == qbpp::equal{0, 2}) は
s が 0 か 2 のときだけ充足されます(qbpp::cons(s, qbpp::equal{0, 2})
と書いても同じです).許容値は任意個・任意の整数を指定できます.
auto e = qbpp::var("e", 5);
// 各頂点に接続する辺を 0 本か 2 本だけ選ぶ
auto deg = 100 * qbpp::cons(qbpp::sum(e) == qbpp::equal{0, 2});
これはグラフの path や cycle を構成する辺を選ぶ問題(各頂点の次数が
0 か 2 のとき充足)などに便利です.許容値がとびとびのため,両側範囲
l <= f <= u では表現できません.制約リストには == {0, 2} と
表示されます.この制約は EasySolver・ExhaustiveSolver・ABS3Solver
で使えます(MIP ソルバーは非対応).double 係数フロントエンド
(DOUBLE_TYPE*)では使えません.
非線形の制約本体
非線形(2 次以上)の式を qbpp::cons() に入れた場合も,等式(x*y + z == 1
など)・範囲(1 <= x*y + z*w <= 2 など)ともにそのまま制約として扱われ,
バンドルされたソルバー(EasySolver / ExhaustiveSolver / ABS3Solver)が
制約を満たすよう探索します.外部の MIP/ILP ソルバーは非線形の制約本体を
受け付けないため,その場合は expand_cons()(後述)で従来のペナルティ式に
展開してから渡してください.
制約が返す値
qbpp::cons() で宣言した制約は,変数割当のもとで違反量の 2 乗を
値として持ちます.制約本体の式を \(f\) とすると,等式制約の値は
範囲制約の値は
です(片側制約は該当する側だけが働きます).重み \(P\) を掛けた制約 \(P \cdot \operatorname{cons}(\cdots)\) はこの値の \(P\) 倍になり, モデル全体の値は
で,バンドルソルバーが報告する Energy と一致します. 全ての制約を満たす解では Energy = objective です.
式の演算規則
制約付きの式 f はモデルの完全な記述です.
f(sol)はソルバーが報告する Energy と一致します.f.cons(sol)は違反している制約の本数を返します(0 なら全充足).目的関数の調整(
+,-, 定数加算)と 0 以外のスカラー倍(重みの一括 スケール),simplify_as_binary(),qbpp::replace()は制約を保ったまま 使えます.f.simplify_as_binary()は目的関数と制約の両方に適用されます. ソルバーに渡す前に 1 回呼んでください — 特にqbpp::replace()で 変数を置換した後に必要です.重みは通常は正の値を使いますが,負の重みも指定できます(制約の減算・ 負号も同様に重みの符号反転として扱われます).負の重みは違反を 「優遇」する特殊な用途向けで,バンドルソルバー (
EasySolver・ExhaustiveSolver・ABS3Solver)でのみ使えます — MIP ソルバー(ハード制約扱い)に渡すとエラーになります.制約の宣言を壊す演算 —
sqr(),式同士の乗算,0 倍(制約が黙って 消えるため),reduce()など — は明示的にエラーになります.
ソルバーごとの意味論
すべてのソルバーが同じ式 f を 1 引数で受け付けます.
ソルバー |
意味論 |
|---|---|
|
ソフト: 制約違反には重みに応じたペナルティが加算され,制約を満たす良い解を効率よく探索する |
|
ソフト: |
|
ハード: 制約は MIP の線形制約として渡される(重みは無視.負の重みの制約が含まれる場合はエラー) |
同一のモデル定義を厳密ソルバーで検証してからヒューリスティックソルバーで スケールアップできます.
#include <qbpp/easy_solver.hpp>
#include <qbpp/exhaustive_solver.hpp>
#include <qbpp/qbpp.hpp>
int main() {
auto z = qbpp::var("z", 4);
qbpp::Expr obj = -3 * z(0) - 5 * z(1) - 2 * z(2) - 7 * z(3);
auto f = obj + 100 * qbpp::cons(z(0) + z(1) + z(2) + z(3) == 2);
f.simplify_as_binary();
qbpp::ExhaustiveSolver exact(f); // ground truth
std::cout << "exact optimum = " << exact.search().energy() << std::endl;
qbpp::EasySolver heuristic(f); // same model, scales up
auto sol = heuristic.search({{"time_limit", 1}});
std::cout << "heuristic = " << sol.energy() << std::endl;
}
ネイティブ制約がある場合,target_energy は「エネルギーが target に達し,
かつ全制約が充足」のときだけ探索を停止します.
EasySolver のデフォルトコールバックはエネルギーと並べて充足の進捗を
表示します.Energy はペナルティ込みの合計,Obj は目的関数部分,
Viol = k/m は m 本の制約のうち k 本が違反中であることを示します.
全制約が充足されると Energy と Obj は一致します.
解の検証
violations() は解に対して全制約を評価し,制約値・境界・違反距離・重みを
報告します.
auto s = qbpp::var("s", 3);
qbpp::Expr obj2 = -1 * s(0) - 2 * s(1) - 3 * s(2);
auto f3 = obj2 + 10 * qbpp::cons(s(0) + s(1) + s(2) == 1);
f3.simplify_as_binary();
qbpp::ExhaustiveSolver solver3(f3);
auto sol3 = solver3.search();
for (const auto& t : f3.violations(sol3)) {
std::cout << t.lower << " <= " << t.value << " <= " << t.upper
<< " violation = " << t.violation << std::endl;
}
std::cout << (f3.is_feasible(sol3) ? "feasible" : "infeasible")
<< std::endl;
従来のペナルティ式への展開
qbpp::expand_cons(f) は,宣言された制約を従来のペナルティ式
(比較演算子で書いた場合と同じ形)に展開した通常の式を
返します.ネイティブ制約に対応しない外部の QUBO/HUBO ツールに渡す場合
などに使います.f 自身を上書きする f.expand_cons() もあります.
展開結果は簡約されていないので,ソルバーに渡す前に
simplify_as_binary() を呼んでください.
auto n0 = qbpp::var("n0"), n1 = qbpp::var("n1"), n2 = qbpp::var("n2");
auto fe = n0 + 10 * qbpp::cons(n1 + n2 == 1);
auto ge = qbpp::expand_cons(fe); // fe is unchanged; ge is a plain Expr
ge.simplify_as_binary();
自由記述ペナルティ
充足のときちょうど値が 0 になる式であれば,従来の QUBO ペナルティ
スタイルの式をそのまま qbpp::cons() に混ぜられます.
auto d = qbpp::var("d"), e = qbpp::var("e"), f = qbpp::var("f");
auto g = qbpp::var("g"), h = qbpp::var("h"), i = qbpp::var("i");
auto mixed = qbpp::cons(100 * (d + e + f == 2) // tracked comparison
+ 200 * (~g * ~h * ~i)); // penalty: not all 0
比較で書いた制約は 1 本ずつ追跡されます.自由記述部分は,その値が 0 の
ときだけ充足と見なされます.デフォルトコールバックはこの部分を
Pen = ...(0 なら充足)として表示し,violations() は境界 [0, 0] の
最終エントリとして報告します.式が非負で最小値 0 になることの保証は
利用者の責任です.
PyQBPP では,式の中の制約部分を qbpp.cons() で囲むと,その部分は
制約とみなされて特別に処理されます.Hi-QUBO にバンドルされている
ソルバーは,宣言された制約を満たすように効率よく探索を行います.
整数線形計画法を cons() で解く
範囲制約と整数線形計画法の求解では,次の整数線形計画問題を, 範囲制約を重み付きのペナルティ式として目的関数に加える方法で解きました:
同じ問題は,制約を qbpp.cons() で作成すると次のように書けます:
import pyqbpp as qbpp
x = qbpp.var("x", between=(0, 10))
y = qbpp.var("y", between=(0, 10))
f = 5 * x + 4 * y
g = -f + 100 * (qbpp.cons(2 * x + 3 * y, between=(0, 24)) +
qbpp.cons(7 * x + 5 * y, between=(0, 54)))
g.simplify_as_binary()
solver = qbpp.EasySolver(g)
sol = solver.search(time_limit=1.0)
print(f"x = {sol(x)}, y = {sol(y)}")
print(f"f = {sol(f)}")
print(f"violated constraints = {g.cons(sol)}")
変更点は,constrain() で作っていたペナルティ式を qbpp.cons() に
置き換えただけです(引数の書き方は同じです).
これだけで 2 つの範囲制約は単なるペナルティ式ではなく制約として宣言され,
ソルバーは制約を満たす解を効率よく探索します.
g.cons(sol) は解 sol で違反している制約の本数を返します(0 なら全制約を充足).
プログラムの出力は以下の通りです:
x = 4, y = 5
f = 40
violated constraints = 0
ナップサック問題の例
もう1つの例として,次のプログラムは簡単なナップサック問題
(容量制約と等式制約)を cons() で解きます:
import pyqbpp as qbpp
x = qbpp.var("x", shape=(6,))
value = [3, 5, 2, 7, 4, 6]
weight = [2, 4, 1, 5, 3, 4]
obj = qbpp.Expr(0)
load = qbpp.Expr(0)
for i in range(6):
obj += -value[i] * x[i]
load += weight[i] * x[i]
f = obj + 100 * qbpp.cons(load, between=(None, 8)) # capacity
f += 10 * qbpp.cons((x[0] + x[1]) == 1) # equality
f.simplify_as_binary()
solver = qbpp.EasySolver(f)
sol = solver.search(time_limit=1)
print("objective =", sol.energy)
従来のペナルティ式からの移行は,制約部分を qbpp.cons() で囲むだけです —
obj + 1000 * (rows + cols) を obj + 1000 * qbpp.cons(rows + cols) に
書き換えます.多くの問題で,同じ制約をペナルティ式のまま解くより
大幅に良い解が得られます.
制約の書き方
制約は qbpp.cons(式 == 整数),または範囲を kwargs で直接指定する
qbpp.cons(式, between=(下限, 上限))(片側は None)で書きます.
重みは制約へのスカラー係数として書き,+ で目的関数や他の制約と
自由に足し合わせられます.
import pyqbpp as qbpp
a, b, c = qbpp.var("a"), qbpp.var("b"), qbpp.var("c")
u, v, w = qbpp.var("u"), qbpp.var("v"), qbpp.var("w")
cons = 1000 * qbpp.cons((a + b + c) == 2) # equality
cons += 1000 * qbpp.cons(10*u + 30*v - 10*w, between=(None, 35)) # one-sided
cons += 500 * qbpp.cons(a + b - c, between=(0, 1)) # two-sided
配列の比較を囲むと要素ごとに 1 本の制約になるので,行列の one-hot 行は 1 文で書けます.
import pyqbpp as qbpp
y = qbpp.var("y", shape=(4, 4))
one_hot = 1000 * qbpp.cons(qbpp.vector_sum(y) == 1) # one per row
制約を蓄積した式に *= を使うと,蓄積済みの全制約の重みを一括で
スケールできます.
import pyqbpp as qbpp
k0, k1, k2 = qbpp.var("k0"), qbpp.var("k1"), qbpp.var("k2")
cons4 = qbpp.cons((k0 + k1 + k2) == 2)
cons4 += qbpp.cons(k0 + k1 - k2, between=(None, 1))
cons4 *= 1000 # scale ALL weights at once
式を print すると目的関数の多項式が表示され,f.cons() は宣言された
制約リストの文字列を返します(重みが 1 のときは係数プレフィックスを
省略,片側制約は片側表示).
import pyqbpp as qbpp
m0, m1, m2 = qbpp.var("m0"), qbpp.var("m1"), qbpp.var("m2")
printed = 1000 * qbpp.cons((m0 + m1 + m2) == 2)
printed += 500 * qbpp.cons(m0 + m1 - m2, between=(0, 1))
print(printed.cons())
出力は次のようになります.
1000 * (m0 +m1 +m2 == 2)
500 * (0 <= m0 +m1 -m2 <= 1)
離散許容値集合
式の値がとびとびの許容値のいずれかに一致することを要求する制約は,
equal=[...] で書けます.qbpp.cons(s, equal=[0, 2]) は s が 0 か 2 の
ときだけ充足されます.許容値は任意個・任意の整数を指定できます.
e = qbpp.var("e", 5)
# 各頂点に接続する辺を 0 本か 2 本だけ選ぶ
deg = 100 * qbpp.cons(qbpp.sum(e), equal=[0, 2])
これはグラフの path や cycle を構成する辺を選ぶ問題(各頂点の次数が
0 か 2 のとき充足)などに便利です.許容値がとびとびのため,両側範囲
between=(l, u) では表現できません.制約リストには == {0, 2} と
表示されます.この制約は EasySolver・ExhaustiveSolver・ABS3Solver
で使えます(MIP ソルバーは非対応).任意精度(pyqbpp.cppint)および
double 係数(pyqbpp.d)バリアントでは使えません.
非線形の制約本体
非線形(2 次以上)の式を qbpp.cons() に入れた場合も,等式(x*y + z == 1
など)・範囲(qbpp.cons(x*y + z*w, between=(1, 2)) など)ともにそのまま
制約として扱われ,バンドルされたソルバー(EasySolver / ExhaustiveSolver /
ABS3Solver)が制約を満たすよう探索します.外部の MIP/ILP ソルバーは非線形の
制約本体を受け付けないため,その場合は expand_cons()(後述)で従来の
ペナルティ式に展開してから渡してください.
制約が返す値
qbpp.cons() で宣言した制約は,変数割当のもとで違反量の 2 乗を
値として持ちます.制約本体の式を \(f\) とすると,等式制約の値は
範囲制約の値は
です(片側制約は該当する側だけが働きます).重み \(P\) を掛けた制約 \(P \cdot \operatorname{cons}(\cdots)\) はこの値の \(P\) 倍になり, モデル全体の値は
で,バンドルソルバーが報告する Energy と一致します. 全ての制約を満たす解では Energy = objective です.
式の演算規則
制約付きの式 f はモデルの完全な記述です.
sol(f)はソルバーが報告する Energy と一致します.f.cons(sol)は違反している制約の本数を返します(0 なら全充足).目的関数の調整(
+,-, 定数加算)と 0 以外のスカラー倍(重みの一括 スケール),simplify_as_binary(),qbpp.replace()は制約を保ったまま 使えます.f.simplify_as_binary()は目的関数と制約の両方に適用されます. ソルバーに渡す前に 1 回呼んでください — 特にqbpp.replace()で 変数を置換した後に必要です.重みは通常は正の値を使いますが,負の重みも指定できます(制約の減算・ 負号も同様に重みの符号反転として扱われます).負の重みは違反を 「優遇」する特殊な用途向けで,バンドルソルバー (
EasySolver・ExhaustiveSolver・ABS3Solver)でのみ使えます — MIP ソルバー(ハード制約扱い)に渡すとRuntimeErrorになります.制約の宣言を壊す演算 —
qbpp.sqr(),式同士の乗算,0 倍(制約が黙って 消えるため),qbpp.reduce()など — はRuntimeErrorになります.
ソルバーごとの意味論
バンドルソルバーは同じ式 f を 1 引数で受け付けます.
ソルバー |
意味論 |
|---|---|
|
ソフト: 制約違反には重みに応じたペナルティが加算され,制約を満たす良い解を効率よく探索する |
|
ソフト: |
外部 MIP ソルバー( |
ハード: 制約は MIP の線形制約として渡される(重みは無視.負の重みの制約が含まれる場合は |
制約付きモデルを外部 MIP ソルバーに渡すには ILP モードを使います —
ilp=True を指定してください(例: qbpp.ScipSolver(f, ilp=True).モデルは
線形である必要があります).ilp=True なしで制約付きモデルを渡すと
エラーになります.
同一のモデル定義を厳密ソルバーで検証してからヒューリスティックソルバーで スケールアップできます.
import pyqbpp as qbpp
z = qbpp.var("z", shape=(4,))
obj = -3*z[0] - 5*z[1] - 2*z[2] - 7*z[3] + qbpp.Expr(0)
f = obj + 100 * qbpp.cons((z[0] + z[1] + z[2] + z[3]) == 2)
f.simplify_as_binary()
exact = qbpp.ExhaustiveSolver(f) # ground truth
print("exact optimum =", exact.search().energy)
heuristic = qbpp.EasySolver(f) # same model, scales up
print("heuristic =", heuristic.search(time_limit=1).energy)
ネイティブ制約がある場合,target_energy は「エネルギーが target に達し,
かつ全制約が充足」のときだけ探索を停止します.
EasySolver のデフォルトコールバックはエネルギーと並べて充足の進捗を
表示します.Energy はペナルティ込みの合計,Obj は目的関数部分,
Viol = k/m は m 本の制約のうち k 本が違反中であることを示します.
全制約が充足されると Energy と Obj は一致します.
解の検証
violations(sol) は解に対して全制約を評価し,制約値・境界・違反距離・
重みの dict を制約ごとに返します.
import pyqbpp as qbpp
s = qbpp.var("s", shape=(3,))
obj = -1*s[0] - 2*s[1] - 3*s[2] + qbpp.Expr(0)
f = obj + 10 * qbpp.cons((s[0] + s[1] + s[2]) == 1)
f.simplify_as_binary()
sol = qbpp.ExhaustiveSolver(f).search()
for t in f.violations(sol):
print(t["lower"], "<=", t["value"], "<=", t["upper"],
" violation =", t["violation"])
print("feasible" if f.is_feasible(sol) else "infeasible")
従来のペナルティ式への展開
qbpp.expand_cons(f) は,宣言された制約を従来のペナルティ式
(比較演算子や qbpp.constrain で書いた場合と同じ形)に展開した通常の
式を返します.ネイティブ制約に対応しない外部の QUBO/HUBO ツールに渡す
場合などに使います.f 自身を上書きする f.expand_cons() もあります.
展開結果は簡約されていないので,ソルバーに渡す前に
simplify_as_binary() を呼んでください.
import pyqbpp as qbpp
n0, n1, n2 = qbpp.var("n0"), qbpp.var("n1"), qbpp.var("n2")
fe = n0 + 10 * qbpp.cons((n1 + n2) == 1)
ge = qbpp.expand_cons(fe) # fe is unchanged; ge is a plain Expr
ge.simplify_as_binary()
自由記述ペナルティ
充足のときちょうど値が 0 になる式であれば,従来の QUBO ペナルティ
スタイルの式をそのまま qbpp.cons() に混ぜられます.
import pyqbpp as qbpp
d, e, f2 = qbpp.var("d"), qbpp.var("e"), qbpp.var("f2")
g, h, i = qbpp.var("g"), qbpp.var("h"), qbpp.var("i")
mixed = qbpp.cons(100 * ((d + e + f2) == 2) # tracked comparison
+ 200 * (~g * ~h * ~i)) # penalty: not all of g,h,i are 0
比較で書いた制約は 1 本ずつ追跡されます.自由記述部分は,その値が 0 の
ときだけ充足と見なされます.デフォルトコールバックはこの部分を
Pen = ...(0 なら充足)として表示し,violations() は境界 [0, 0] の
最終エントリとして報告します.式が非負で最小値 0 になることの保証は
利用者の責任です.