モチベーション
「あたらしい数理最適化 python言語とGurobiで解く」
を教科書として、数理最適化の基本となるところを、肌感覚で理解したい。そのために、
- 基礎的な内容
- その例題演習
を書くことで、自身の理解と、新たな仕事の確保を目指します。私の大本のモチベーションはここを参照してください。
この記事を読むとできること
- Google OR-Toolで制約条件付きの簡単な問題を解けるようにする
を目指します。
ステップ1: 写経
例の如く、Python用のStart Guideから、写経をしていきます。以下、引用。
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(-solver.infinity(), solver.infinity(), 'x')
y = solver.NumVar(-solver.infinity(), solver.infinity(), 'y')
# Constraint 1: x + 2y <= 14.
constraint1 = solver.Constraint(-solver.infinity(), 14)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 2)
# Constraint 2: 3x - y >= 0.
constraint2 = solver.Constraint(0, solver.infinity())
constraint2.SetCoefficient(x, 3)
constraint2.SetCoefficient(y, -1)
# Constraint 3: x - y <= 2.
constraint3 = solver.Constraint(-solver.infinity(), 2)
constraint3.SetCoefficient(x, 1)
constraint3.SetCoefficient(y, -1)
objective = solver.Objective()
objective.SetCoefficient(x, 1)
objective.SetCoefficient(y, 1)
objective.SetMaximization()
# Call the solver and display the results.
solver.Solve()
print('Solution:')
print('x = ', x.solution_value())
print('y = ', y.solution_value())
if __name__ == '__main__':
main()
今回も、実行すればちゃんと最後まで通りました。ありがてぇ。
ステップ2: 意味を理解する
さて、写経は済んだので意味を理解していきましょう。
ソルバのインスタンス
solver = pywraplp.Solver(
'SolveSimpleSystem', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
これは明らかにソルバのインスタンスを立てていますね。前回同様、第一引数は、名前です。前回は調べなかった、pywraplp.Solver.GLOP_LINEAR_PROGRAMMINGは何を意味するんでしょうか。GLOPについて調べたところ、こんなページを見つけました。どうやら、Googleが作った線形ソルバのようです。計算量も空間量も優れていると書かれていますが、何がどうすごいのかまでは書いていませんね...。とはいえ、このページの最後にThe Stigler dietなる線形計画問題を見つけました。どうやらこのおっさんが提唱した、人が最低限の栄養を確保した際にかかるコストの最小値を問う問題のようです。
変数
# Create the variables x and y.
x = solver.NumVar(-solver.infinity(), solver.infinity(), 'x')
y = solver.NumVar(-solver.infinity(), solver.infinity(), 'y')
変数の設定。前回と違って、無限の定義域が設定されています。Cでよくある持ち方だと思いますが、solver.infinty()が無限を意味するようですね。
制約条件
# Constraint 1: x + 2y <= 14.
constraint1 = solver.Constraint(-solver.infinity(), 14)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 2)
これが制約条件ですね（Constraintは制約を意味します）。制約条件を作る際に、その条件の値域を設定して、そこに係数付きで変数を放り込むようです。今回は不等式ですが、等式の場合はどうするんでしょうね。
# Constraint 1: x + 2y = 14.
constraint1 = solver.Constraint(14, 14)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 2)
とすると、誤差が引っかかって制約条件を満たさない気がします。試してみましょう。まずはそのままの場合。
Solution:
x = 5.999999999999998
y = 3.9999999999999996
>>>
明らかに誤差がついてきています。これは期待大。上式にしてやってみます。
Solution:
x = 6.0
y = 4.0
>>>
やべぇ、ちゃんと解いた...。等式にすることでむしろ精度があがってます。制約がキツいせいで、逆に大域的最適解に素早く収束したんでしょうか。もうちょっときつくしてみましょうか。
# Constraint 3: x - y <= 2.
constraint3 = solver.Constraint(-solver.infinity(), 2)
constraint3.SetCoefficient(x, 1)
constraint3.SetCoefficient(y, -1)
を、
# Constraint 3: x - y = 2.
constraint3 = solver.Constraint(2, 2)
constraint3.SetCoefficient(x, 1)
constraint3.SetCoefficient(y, -1)
に変えてみましょう。完全に連立方程式を解かせていますが、あんまりきにしない。
Solution:
x = 6.0
y = 4.0
>>>
ふむ。それでこそソルバだという気分にさせられます。よくできていて、面白い。しかし、なんで条件をきつくすると精度も上がるんですかね....。等式という強い制約条件を満たせるから（実行可能解があるから）さらに正確に解いた、が正解でしょうか。
なお、等式の取扱についてgoogle先生に聞いたところ、Constraint Optimizationなるものを教えてくれました。これは目的関数云々より、制約条件を満たす変数の組を探索してくれるもの。実行可能性を探してくれるともいえそうです。
目的関数
objective = solver.Objective()
objective.SetCoefficient(x, 1)
objective.SetCoefficient(y, 1)
objective.SetMaximization()
ここは前と変わらず。大体わかってきました。
最適化の実行
solver.Solve()
print('Solution:')
print('x = ', x.solution_value())
print('y = ', y.solution_value())
ここも疑問点はありませんね。pythonの、値渡しと参照渡しを同時に行う仕様がバリバリ使われていて、気を使うところでもありますね。値渡ししかしない世界に長いこと居すぎたので、下手に書くと偉いことになりそうです。
まとめと考察
とりあえず、制約条件つきの線形計画問題は解くことができました。だいたいの構造も把握してきたので、そろそろ実践問題を一個解いてみたいですね。調べてる途中で行き着いた、以下の話はやって見る価値がありそうです。
- The Stigler diet
- Constraint Optimization
- その他、私が知らない最適化問題の形式
終わりに
自分の記事のクオリティを高めたいので、質問でも、次の要望、ただの感想、何でもいいのでコメントを残して頂けると幸いです。ちゃんと答えますので、よしなしによろしくお願いいたします。