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

整数線形計画法の求解

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

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

この問題の最適解は $x=4$, $y=5$ であり、目的関数値は $40$ です。

以下のQUBO++プログラムは、整数変数 x，y を使い，Easy Solverでこの最適解を求めます:

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

int main() {
  auto x = 0 <= qbpp::int_var("x") <= 10;
  auto y = 0 <= qbpp::int_var("y") <= 10;
  auto f = 5 * x + 4 * y;
  auto g = -f + 100 * qbpp::cons(2 * x + 3 * y <= 24) +
           100 * qbpp::cons(7 * x + 5 * y <= 54);
  g.simplify();
  std::cout << "g = " << g << std::endl;

  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;
}

このQUBO++プログラムでは、

f は目的関数を表し、
qbpp::cons(2 * x + 3 * y <= 24) と qbpp::cons(7 * x + 5 * y <= 54) は2つの制約を表し、
g はこれらを1つの最適化式にまとめたものです。

目標が最大化であるため、目的関数は -f として符号を反転しています。 qbpp::cons() で囲んだ比較式は制約として特別に処理され、バンドルされたソルバーは制約を満たすように効率よく探索を行います。制約に違反すると，違反量（この例では左辺が右辺を超えた量）の 2 乗に重み 100 を掛けたペナルティが加わります．重みは，制約を破って目的関数を改善しても得にならないよう，十分大きな値にします．

g に対してEasy Solverインスタンスを作成し、制限時間1.0秒で探索を実行します。 g.cons(sol) は解 sol で違反している制約の本数を返します（0 なら全制約を充足）。プログラムの出力は以下の通りです:

g = -5*x -4*y
+100*cons(2*x +3*y <= 24)
+100*cons(7*x +5*y <= 54)
x = 4, y = 5
f = 40
violated constraints = 0

g を出力すると，目的関数の多項式に続いて，宣言した制約が 1 行ずつ表示されます．ソルバーが正しく最適解を見つけたことが確認できます。

範囲制約

制約には，片側の不等式のほかに，範囲制約 $l \leq f \leq u$ も書けます． qbpp::cons() の中に書ける比較式は次の 4 種類です:

書き方	意味
`qbpp::cons(f <= u)`	$f \leq u$
`qbpp::cons(f >= l)`	$f \geq l$
`qbpp::cons(l <= f <= u)`	$l \leq f \leq u$（範囲制約）
`qbpp::cons(f == v)`	$f = v$

例えば，上の問題に「$x - y$ が 1 以上 3 以下」という範囲制約を加えると，最適解は $x=5$，$y=3$（目的関数値 $37$）に変わります:

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

int main() {
  auto x = 0 <= qbpp::int_var("x") <= 10;
  auto y = 0 <= qbpp::int_var("y") <= 10;
  auto f = 5 * x + 4 * y;
  auto g = -f + 100 * qbpp::cons(2 * x + 3 * y <= 24) +
           100 * qbpp::cons(7 * x + 5 * y <= 54) +
           100 * qbpp::cons(1 <= x - y <= 3);
  g.simplify();

  auto solver = qbpp::ExhaustiveSolver(g);
  auto sol = solver.search();
  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;
}

ここでは Exhaustive Solver で全探索し，厳密な最適解を求めています．プログラムの出力は以下の通りです:

x = 5, y = 3
f = 37
violated constraints = 0

cons() で宣言した制約は、バンドルされた 3 つのソルバー（Easy Solver・ Exhaustive Solver・ABS3 Solver）すべてで同じ意味を持ち、Exhaustive Solver はこのペナルティ込みエネルギーの厳密な最小解を返します。 MIP ソルバー（SCIP 等）ではハード制約（必ず満たすべき制約）として扱われます。詳細はネイティブ制約をご覧ください。

ペナルティ式による範囲制約

qbpp::cons() で囲まない比較式は，QUBO・HUBO として表せるペナルティ式（制約が満たされるとき，かつそのときに限り最小値 0 をとる式）になります．バイナリ変数しか扱えない外部の QUBO ソルバーに渡すときは，このペナルティ式を使います．

範囲制約 $l\leq f\leq u$ のペナルティ式は，範囲 $[l,u]$ の値をとる補助整数変数 $a$ を導入して作ります:

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

この式 $g$ は $f=a$ のときに限り最小値0をとります。 $a$ は $[l,u]$ の任意の整数値をとれるため、式 $g$ は $f$ 自身が同じ範囲内の整数値をとる場合に限り0を達成します。補助変数 $a$ はバイナリエンコードの整数変数として導入されます． $f$ が線形式の場合、$g$ はQUBO式になります。 $f$ が2次以上の場合、$g$ はHUBO式になります（$g$ の次数は $f$ の2倍になるため）。

ペナルティ式では範囲制約ごとにこの補助変数が必要ですが、 cons() で宣言した場合は補助変数を導入しません。モデルの変数数と項数がどれだけ変わるかは、ネイティブ制約の「ペナルティ式との比較」で比べています。 cons() で書いたモデルは，qbpp::expand_cons() でペナルティ式に展開できます．

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