範囲制約と整数線形計画法の求解¶
範囲制約の多項式定式化
\(f\) をバイナリ変数の多項式とします. 範囲制約は \(l<u\) に対して \(l\leq f\leq u\) の形式を持ちます. 目標は,範囲制約が満たされる場合に限り最小値0をとる多項式を設計することです.
鍵となるアイデアは,範囲 \([l,u]\) の値をとる補助整数変数 \(a\) を導入することです. 以下の式を考えます:
この式 \(g\) は \(f=a\) のときに限り最小値0をとります. \(a\) は \([l,u]\) の任意の整数値をとれるため,式 \(g\) は \(f\) 自身が同じ範囲内の整数値をとる場合に限り0を達成します.
この補助変数の手法を用いて,Hi-QUBOは範囲制約を実装しています. \(f\) が線形式の場合,\(g\) はQUBO式になります. \(f\) が2次以上の場合,\(g\) はHUBO式になります（\(g\) の次数は \(f\) の2倍になるため）.
NOTE Hi-QUBOは内部的に軽量な改善を施しており,範囲制約をわずかに少ないバイナリ変数数で符号化できます. 詳細は比較演算子に記載されています.
整数線形計画法の求解
整数線形計画法のインスタンスは,目的関数と複数の線形制約から構成されます. 例えば,以下の整数線形計画問題は2つの変数,1つの目的関数,2つの制約を持ちます:
この問題の最適解は \(x=4\), \(y=5\) であり,目的関数値は \(40\) です.
以下のHi-QUBOプログラムは,Easy Solverを使用してこの最適解を求めます:
#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 * (c1 + 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 << "c1 = " << sol(c1) << ", c2 = " << sol(c2) << std::endl;
std::cout << "c1.body(sol) = " << c1.body(sol) << ", c2.body(sol) = " << c2.body(sol) << std::endl;
}
このHi-QUBOプログラムでは,
fは目的関数を表し,c1とc2は範囲制約を表し,gはこれらを1つの最適化式にまとめたものです.
目標が最大化であるため,目的関数は -f として符号を反転しています.
制約 c1 と c2 には重み100のペナルティを付け,高い優先度で制約が満たされるようにしています.
g に対してEasy Solverインスタンスを作成し,制限時間1.0秒で探索を実行します.
最適解 sol を取得した後,x,y,f,c1,c2,c1.body(sol),c2.body(sol) の値を出力します.
プログラムの出力は以下の通りです:
x = 4, y = 5
f = 40
c1 = 0, c2 = 0
c1.body(sol) = 23, c2.body(sol) = 53
ここで,
c1は制約0 <= 2 x + 3 y <= 24のペナルティ式（制約が満たされると 0）であり,c1.body()は線形式2 x + 3 yを返し,c1.body(sol)はその値をsolで評価します.
ソルバーが正しく最適解を見つけたことが確認できます.
ネイティブ制約 cons() による記述
上のプログラムでは,制約 c1 と c2 を重み付きのペナルティ式として目的関数に加えました.
Hi-QUBO では,さらに一歩進めて,制約を 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)) に書き換えただけです.
qbpp::cons() で囲まれた部分は制約として特別に処理され,
バンドルされたソルバーは宣言された制約を満たすように効率よく探索を行います.
g.cons(sol) は解 sol で違反している制約の本数を返します（0 なら全制約を充足）.
プログラムの出力は以下の通りです:
x = 4, y = 5
f = 40
violated constraints = 0
cons() で宣言した制約は,バンドルされた 3 つのソルバー（Easy Solver・
Exhaustive Solver・ABS3 Solver）すべてで同じ意味を持ち,Exhaustive Solver は
このペナルティ込みエネルギーの厳密な最小解を返します.
MIP ソルバー（Gurobi 等）ではハード制約（必ず満たすべき制約）として扱われます.
詳細はネイティブ制約をご覧ください.
範囲制約の多項式定式化
\(f\) をバイナリ変数の多項式とします. 範囲制約は \(l<u\) のもとで \(l\leq f\leq u\) の形式を持ちます. 目標は,範囲制約が満たされるときかつそのときに限り最小値0を取る多項式を設計することです.
鍵となるアイデアは,範囲 \([l,u]\) の値を取る補助整数変数 \(a\) を導入することです. 以下の式を考えます：
この式 \(g\) は \(f=a\) のときちょうど最小値0を取ります. \(a\) は \([l,u]\) の任意の整数値を取れるため, 式 \(g\) が0になるのは \(f\) 自体が同じ範囲内の整数値を取るときかつそのときに限ります.
この補助変数の手法を用いて,PyQBPPは constrain() 関数により範囲制約を実装しています.
NOTE PyQBPPは内部的に軽量な改善を施しており,範囲制約をわずかに少ないバイナリ変数数で符号化できます. 詳細は比較演算子に記載されています.
整数線形計画法の求解
整数線形計画法のインスタンスは,目的関数と複数の線形制約から構成されます. 例えば,以下の整数線形計画は2つの変数,1つの目的関数,2つの制約を持ちます：
この問題の最適解は \(x=4\), \(y=5\) で,目的関数の値は \(40\) です.
以下のPyQBPPプログラムは,Easy Solverを使ってこの最適解を求めます：
import pyqbpp as qbpp
x = qbpp.var("x", between=(0, 10))
y = qbpp.var("y", between=(0, 10))
f = 5 * x + 4 * y
c1 = qbpp.constrain(2 * x + 3 * y, between=(0, 24))
c2 = qbpp.constrain(7 * x + 5 * y, between=(0, 54))
g = -f + 100 * (c1 + c2)
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"c1 = {sol(c1)}, c2 = {sol(c2)}")
print(f"2x+3y = {sol(c1.body)}, 7x+5y = {sol(c2.body)}")
このプログラムでは,
fは目的関数を表し,c1とc2はconstrain()を使って作成された範囲制約を表し,gはそれらを1つの最適化式にまとめたものです.
目標が最大化であるため,目的関数は -f として符号を反転しています.
制約 c1 と c2 は重み100のペナルティを付けて,高い優先度で満たされるようにしています.
g に対してEasy Solverのインスタンスを作成し,制限時間1.0秒を search() のパラメータとして渡して探索を実行します.
最適解 sol を得た後,プログラムは x,y,f,c1,c2,および制約本体の式の値を出力します.
プログラムの出力は以下の通りです：
x = 4, y = 5
f = 40
c1 = 0, c2 = 0
c1.body(sol) = 23, c2.body(sol) = 53
ここで,
c1は制約0 <= 2x + 3y <= 24のペナルティであり,c1.bodyは線形式2x + 3yを表します.
ソルバーが正しく最適解を見つけていることが確認できます.
ネイティブ制約 cons() による記述
上のプログラムでは,制約 c1 と c2 を重み付きのペナルティ式として目的関数に加えました.
PyQBPP では,さらに一歩進めて,制約を 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() を使うことだけです（引数の書き方は同じで,
範囲制約は between=(l, u),等式制約は equal=n で指定します）.
qbpp.cons() で作成した部分は制約として特別に処理され,
バンドルされたソルバーは宣言された制約を満たすように効率よく探索を行います.
g.cons(sol) は解 sol で違反している制約の本数を返します（0 なら全制約を充足）.
プログラムの出力は以下の通りです:
x = 4, y = 5
f = 40
violated constraints = 0
cons() で宣言した制約は,バンドルされた 3 つのソルバー（Easy Solver・
Exhaustive Solver・ABS3 Solver）すべてで同じ意味を持ち,Exhaustive Solver は
このペナルティ込みエネルギーの厳密な最小解を返します.
MIP ソルバー（Gurobi 等）では ilp=True を指定するとハード制約
（必ず満たすべき制約）として扱われます.
詳細はネイティブ制約をご覧ください.