範囲制約と整数線形計画法の求解

範囲制約の多項式定式化

\(f\) をバイナリ変数の多項式とします. 範囲制約は \(l<u\) に対して \(l\leq f\leq u\) の形式を持ちます. 目標は,範囲制約が満たされる場合に限り最小値0をとる多項式を設計することです.

鍵となるアイデアは,範囲 \([l,u]\) の値をとる補助整数変数 \(a\) を導入することです. 以下の式を考えます:

\[ \begin{aligned} g &= (f-a)^2 \end{aligned} \]

この式 \(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つの制約を持ちます:

\[\begin{split} \begin{aligned} \text{Maximize: } & & & 5x + 4y \\ \text{Subject to: } & && 2x + 3y \le 24 \\ & & & 7x + 5y \le 54 \end{aligned} \end{split}\]

この問題の最適解は \(x=4\), \(y=5\) であり,目的関数値は \(40\) です.

以下のHi-QUBOプログラムは,Easy Solverを使用してこの最適解を求めます:

range-constraints-and-ilp-program1.cpp
#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() で囲んで制約であることを明示できます:

range-constraints-and-ilp-program2.cpp
#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\) を導入することです. 以下の式を考えます：

\[ \begin{aligned} g &= (f-a)^2 \end{aligned} \]

この式 \(g\) は \(f=a\) のときちょうど最小値0を取ります. \(a\) は \([l,u]\) の任意の整数値を取れるため, 式 \(g\) が0になるのは \(f\) 自体が同じ範囲内の整数値を取るときかつそのときに限ります.

この補助変数の手法を用いて,PyQBPPは constrain() 関数により範囲制約を実装しています.

NOTE PyQBPPは内部的に軽量な改善を施しており,範囲制約をわずかに少ないバイナリ変数数で符号化できます. 詳細は比較演算子に記載されています.

整数線形計画法の求解

整数線形計画法のインスタンスは,目的関数と複数の線形制約から構成されます. 例えば,以下の整数線形計画は2つの変数,1つの目的関数,2つの制約を持ちます：

\[\begin{split} \begin{aligned} \text{Maximize: } & & & 5x + 4y \\ \text{Subject to: } & && 2x + 3y \le 24 \\ & & & 7x + 5y \le 54 \end{aligned} \end{split}\]

この問題の最適解は \(x=4\), \(y=5\) で,目的関数の値は \(40\) です.

以下のPyQBPPプログラムは,Easy Solverを使ってこの最適解を求めます：

range-constraints-and-ilp-program1.py
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() で作成して制約であることを明示できます:

range-constraints-and-ilp-program2.py
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 を指定するとハード制約 （必ず満たすべき制約）として扱われます. 詳細はネイティブ制約をご覧ください.