この記事を読むとできること
- Google OR-Toolsを使って整数計画問題が解けるようにする
を目指します。参考図書として以下を利用します。
「あたらしい数理最適化 python言語とGurobiで解く」
イントロ
前回まで取り扱っていた問題は線形計画問題でした。今回は線形ではなく、整数計画問題について取り組みたいと思います。
整数計画問題とは
今まで取り扱っていた線形計画問題は、例えば、以下のような定式化がされていました。
\begin{align}
&\rm{maximize} & & y+z \\
&\rm{subject\ to} & & x+2y+z \leq 30\\
& & & x+4y+8z \leq 80\\
& & & x, y, z\geq0
\end{align}
この場合、3つの変数$x, y, z$は実数であればなんでもよかったのが、線形計画問題でした。ちょっと解いてみましょう。
from __future__ import print_function
from ortools.linear_solver import pywraplp
def main():
solver = pywraplp.Solver(
'SolveSimpleSystem', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
# Create the variables x and y.
x = solver.NumVar(0, solver.infinity(), 'x')
y = solver.NumVar(0, solver.infinity(), 'y')
z = solver.NumVar(0, solver.infinity(), 'z')
# Constraint 1: x + 2y + z <= 30.
constraint1 = solver.Constraint(0, 30)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 2)
constraint1.SetCoefficient(z, 1)
# Constraint 2: x + 4y +8z <= 80
constraint2 = solver.Constraint(0, 80)
constraint2.SetCoefficient(x, 1)
constraint2.SetCoefficient(y, 4)
constraint2.SetCoefficient(z, 8)
# objective: y+z
objective = solver.Objective()
objective.SetCoefficient(y, 1)
objective.SetCoefficient(z, 1)
objective.SetMaximization()
# Call the solver and display the results.
solver.Solve()
print('Solution:')
print('x = ', x.solution_value())
print('y = ', y.solution_value())
print('z = ', z.solution_value())
if __name__ == '__main__':
main()
実行結果はこちらです。
x = 0.0
y = 13.333333333333332
z = 3.333333333333335
どうやら$x=0, y=13\frac{1}{3},z=3\frac{1}{3}$が答えのようですね。今回取り扱う整数計画問題は、以下のように変数が整数でなければならないという制約があるものです。
\begin{align}
&\rm{maximize} & & y+z \\
&\rm{subject\ to} & & x+2y+z \leq 30\\
& & & x+4y+8z \leq 80\\
& & & x, y,z は非負の整数
\end{align}
目的関数がどうとか、制約条件に整数しか使ってはいけないとか言うわけではありません。あくまで変数が整数である、という制約を持つということです。線形計画問題としては、$y, z$はともに実数でしたから、整数ではありません。
ガイドに沿ってやってみる
さて、整数計画問題についてもガイドが書かれています。googleの中の人は優秀だぁ...。これを参考に、まずは使い方を学習していきたいと思います。
本当に実行可能解があるのか
さて、整数計画問題は、線形計画問題に対して、さらなる制約を加えたものだと述べました。("整数"と言う制約に限りませんが)制約が加わっているので、実際に制約を満たす変数の組み合わせ（＝実行可能解）があるとは限りません。そこで、それを知るために、Constraint programming (CP) solverが存在します。このページで詳細が書かれているので、やってみましょう。(以下、引用)
from __future__ import print_function
from ortools.sat.python import cp_model
def main():
model = cp_model.CpModel()
var_upper_bound = max(50, 45, 37)
x = model.NewIntVar(0, var_upper_bound, 'x')
y = model.NewIntVar(0, var_upper_bound, 'y')
z = model.NewIntVar(0, var_upper_bound, 'z')
model.Add(2*x + 7*y + 3*z <= 50)
model.Add(3*x - 5*y + 7*z <= 45)
model.Add(5*x + 2*y - 6*z <= 37)
model.Maximize(2*x + 2*y + 3*z)
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
print('Maximum of objective function: %i' % solver.ObjectiveValue())
print()
print('x value: ', solver.Value(x))
print('y value: ', solver.Value(y))
print('z value: ', solver.Value(z))
if __name__ == '__main__':
main()
前回と色々違いますが、model.NewIntVarを指定することで整数計画問題を指定しているようですね。また、model.addで直接、制約条件を挿入できるようですね。全体の構造として、modelを作り、それをsolverに食わせることで、解いているようです。実際に実行してみましょう。
Maximum of objective function: 35
x value: 7
y value: 3
z value: 5
ちゃんと、整数計画問題が解けました！それでは、今の問題に当てはめて解いてみましょう。今回の問題はこれでしたね。
\begin{align}
&\rm{maximize} & & y+z \\
&\rm{subject\ to} & & x+2y+z \leq 30\\
& & & x+4y+8z \leq 80\\
& & & x, y,z は非負の整数
\end{align}
早速コードに落としていきます。 以下のようになりました。
これを実行してみましょう！
from __future__ import print_function
from ortools.sat.python import cp_model
def main():
model = cp_model.CpModel()
var_upper_bound = 100
x = model.NewIntVar(0, var_upper_bound, 'x')
y = model.NewIntVar(0, var_upper_bound, 'y')
z = model.NewIntVar(0, var_upper_bound, 'z')
model.Add(1 * x + 2 * y + 1 * z <= 30)
model.Add(1 * x - 4 * y + 8 * z <= 80)
model.Maximize(y + z)
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
print('Maximum of objective function: %i' % solver.ObjectiveValue())
print()
print('x value: ', solver.Value(x))
print('y value: ', solver.Value(y))
print('z value: ', solver.Value(z))
if __name__ == '__main__':
main()
その結果、以下のように解くことができました。
Maximum of objective function: 22
x value: 0
y value: 8
z value: 14
注意点
実行可能解の集合が凸型のもの（すべての制約条件を満たさなければならないタイプはこれに該当します。）に対して、MIPソルバ（後半で説明します）のほうが高速です。ただし、実行可能解の集合が、凸型でない場合（例えば、3つある制約条件のうち、1つを満たせば良い、といったorで結ばれた制約条件を抱えている場合）は、CPソルバのほうが高速です。実世界の問題だと、凸型でなくて当たり前なんじゃないかな...。
振り返り
今回は、整数計画問題を解くために、CPソルバを利用しました。この記事の上の方でも紹介しましたが、MIPソルバも利用することがあります。MIPソルバを説明するために、少し話をしましょう。最適化の世界には、整数計画問題と線形計画問題の拡張として、混合整数計画問題が存在します。先程出てきたMIPソルバはそちらを解くためのソルバで、もちろん整数計画問題も解けます。MIPソルバは次回あたり、混合整数計画問題は次次回あたりに述べることとします。
終わりに
自分の記事のクオリティを高めたいので、質問でも、次の要望、ただの感想、何でもいいのでコメントを残して頂けると幸いです。