(混合)整数計画問題ソルバーで遊ぶ #kuac2021

この記事は「Kyoto University Advent Calendar 2021」の19日目の記事です。
adventar.org

はじめに

こんにちは．

ほとんどの皆さん初めまして．電タクと言います．工学部電気電子工学科の4回生です．

今回はPython-MIPを使って整数計画問題を解いていきましょう．

本稿は読者の予備知識はそこまで仮定しません．連立方程式が解けて何かしらプログラムを書いたことがあれば理解できると思います．

はじめに
- もくじ
整数計画問題とは
人はなぜ整数計画問題を解くのか？それはそこに整数計画問題があったから
Python-MIPを使ってみる
解きたい問題の説明
(0-1)整数計画問題に定式化してみる
あとはソルバーが勝手に解いてくれる
さいごに

整数計画問題とは

まず初めに整数計画問題について説明します．詳しい説明は計算困難問題に対するアルゴリズム理論や組合せ最適化などの良著に譲ることにして，ここでは簡単な説明にとどめます．

整数計画問題(IP; Integer Programming)と言うのは，制約を満たしながら変数が整数値をとるように動かして目的の関数の値を最大化，あるいは最小化する問題です．

これを数式を使って言い直すとこんな感じです．変数を x1，x2，...，xn とします．x1，x2，...，xn は整数値のみを取ります．目的関数を f とします．f は整数の定数 a1，a2，...，an を用いて f = a1 * x1 + a2 * x2 + ... + an * xn = Σ_i ai * xi (1 ≤ i ≤ n)とかけます．次に制約を定めます．m × n 個の整数の定数 b11，b12，...，b1n，b21，b22，...，b2n，...，bm1，bm2，...，bmn と m 個の整数の定数 c1，c2，...，cm を用いて Σ_i bij * xi ≤ cj (1 ≤ i ≤ n，1 ≤ j ≤ m) を必ず満たすこととします．これらの m 個の不等式制約のもとで目的関数 f を最大化，または最小化します．

とてもわかりにくいですね．具体例で見てみましょう．

(TODO: 後でかく)
(追記：どうせ下で実際に定式化して解くのでそっち見てもらいたいです．疲れが取れたら書くかも．ごめんなさい．)

人はなぜ整数計画問題を解くのか？それはそこに整数計画問題があったから

さて，次になぜ整数計画問題を解くことができると嬉しいのかという話を簡単にします．

簡単に言うと整数計画問題はNP困難であり，また高性能なソルバーを利用して簡単に厳密解を求めることができるのが整数計画問題を解く利点です．つまりどういうことかというのを以下で説明します．

まずNP困難というのはどういうことかですね．NP困難とは計算理論の分野で問題をクラス分けした時の問題のクラスの一つです．長くなるので厳密さを排していうと，NP困難な問題が高速に解けると大体の問題がその解き方を利用して解くことができるようになります．まあ細かい話は計算理論の基礎の3巻とかを読んでいただければ良いです．大事なことは，(計算に時間がかかるものを含む)大体の問題がNP困難問題を解くことで同時に解くことができるということです．つまり，なにか解くのが難しそうな問題でも変形してNP困難な問題にしてしまえばNP困難問題のソルバーに代わりに解いてもらうことで元の問題の答えをゲットできるということです．余談ですが，NP困難な問題を(多項式時間で)高速に解くアルゴリズムは今現在見つかっていません．もし見つけたという方は後でこっそり教えてください．

次にソルバーの話ですね．ソルバーとはその問題の入力を与えると解いて答えを出力してくれるプログラムのことです．仰々しい名前ですが，普通に問題を解くプログラムはソルバーだと思ってもいいと思います．

さて，解きたい問題をNP困難な問題に変形したとして，結局そのNP困難な問題のソルバーがなければ答えを得ることができません．ですが，整数計画問題では優秀なソルバーを簡単に利用することができます．有名なソルバーにCPLEXやGurobi，COIN-CBCなどがあります．またそれらのソルバーを簡単に利用できるようにPythonからそれらを使えるインターフェースとしてPuLPやPython-MIPなどがあります．今回はPython-MIPを利用して解いてみようと思います．

Python-MIPを使ってみる

まずPython-MIPを利用するための準備をしましょう．自分の環境はMacですが大体同じ感じで行けるはずです．

電タクの環境

macOS Monterey version 12.0.1

Python 3.9.7

pip 21.3.1

まずPythonはありますか？なければ入れましょう．terminal.appで $ 以降のコマンドを叩いてみて Python 数字.数字.数字って表示されればokです．

$ python3 --version

Python 3.9.7

次に，pipはインストールされていますか？terminal.appで $ 以降のコマンドを叩いてみて何かしらバージョンが表示されればokです．

$ pip --version

pip 21.3.1 from /usr/hoge/pip (python 3.9)

なければ公式のやり方にしたがって入れましょう．

$ curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py
$ python3 get-pip.py

インストール後バージョンが表示されるようになったらokです．
最後にPython-MIPをpipでインストールしましょう．公式の通り，以下のコマンドを叩くだけです．

$ pip install mip

さて，これでPython-MIPを解く準備は出来ました．では何か整数計画問題を解いてみましょう．

解きたい問題の説明

今回は次のようなゲームでどんな装備が最強かというのを考えてみましょう．(電タクは最近モンハンにはまっています)

プレイヤーは頭・体に一つずつ装備を付けられます
プレイヤーは武器を一つ装備できます
各装備，武器には物理・魔法・速度・防御力・ランクそれぞれの値が設定されています．
プレイヤーはレベルが設定されていて，各装備はランクの合計がレベル以下になるように装備しないといけません

こういうゲームがあった時に，どんな装備を選べば最強かを考えましょう．

ゲームでは装備を全て紙に書き出してあれやこれや考えるのも楽しいですが，効率重視派の人はシミュレータを使って装備を決める人もいると思います．つまりシミュレータの実装をしてみようという話ですね．

まあこれくらいだと頭，体，武器の全てを全探索しても多分間に合うのですが，アップデートで今後さらに装備やルールが複雑になるかも知れません．そうなった時を見越して整数計画問題として解いておけばアップデートのたびに一から書き直す必要がなくなるわけですね．

(0-1)整数計画問題に定式化してみる

まずこの問題で決めないといけないのは，どの装備を選ぶかということですね．ですので次のような l + m + n 個の変数 a1，a2，....，al，b1，b2，....，bm，c1，c2，....，cn を用意します．

l 個の頭装備のうち．ai = 1 の時 i番目の頭装備を装備し ai = 0 の時 i番目の頭装備を装備しない．ai は0以上1以下の整数値をとる．
m 個の体装備のうち．bi = 1 の時 i番目の体装備を装備し bi = 0 の時 i番目の体装備を装備しない．bi は0以上1以下の整数値をとる．
n 個の武器のうち．ci = 1 の時 i番目の武器を装備し ci = 0 の時 i番目の武器を装備しない．ci は0以上1以下の整数値をとる．

これらの変数が満たすべき制約を考えていきましょう．以下ルールを前章で定めた数字で呼びます．

ルール1より，頭装備と体装備は一つずつしか装備できません．これは以下の制約で表現できます．

$\sum_{i=1}^l ai \leq 1$
$\sum_{i=1}^m bi \leq 1$

これは簡単ですね．一つでも装備すれば左辺は1になるので，二つ以上は装備できないことになります．

同様に武器を表す変数 ci も以下の制約を満たします．

$\sum_{i=1}^n ci \leq 1$

次にルール3を受けて，自分が装備した装備の物理，魔法，速度，防御力，ランクの合計を表す変数を導入しましょう．それぞれ P，M，S，D，R としましょう．装備 a，b，cの i 番目のステータスを物理 pai，pbi，pci，魔法 mai，mbi，mci，.... のように定めると(伝われ)，装備の合計は次のように計算することができます．(mが重複していてわかりづらいです，ごめんなさい💦)

$P = \sum_{i=1}^l ai * pai + \sum_{i=1}^m bi * pbi + \sum_{i=1}^n ci * pci$
$M = \sum_{i=1}^l ai * mai + \sum_{i=1}^m bi * mbi + \sum_{i=1}^n ci * mci$
$S = \sum_{i=1}^l ai * sai + \sum_{i=1}^m bi * sbi + \sum_{i=1}^n ci * sci$
$D = \sum_{i=1}^l ai * dai + \sum_{i=1}^m bi * dbi + \sum_{i=1}^n ci * dci$
$R = \sum_{i=1}^l ai * rai + \sum_{i=1}^m bi * rbi + \sum_{i=1}^n ci * rci$

これは装備しているもののみ変数が 1 をとるためこの式によって計算できます．

さいごにルール4の制約を考えましょう．プレイヤーのレベルをLebel
とすると，以下の制約が必要になります．

$R \leq Lebel$

これはそのままですね．

さいごに目的関数を考えましょう．最強装備を求めるということでしたが，装備の組み合わせによって変動するステータスが多くて簡単には順序づけすることができません．今回は割と現実的な目的関数として，物理，魔法，速度の最低ラインを定めた時の防御力最強装備を求めましょう．

ユーザが自分で決めた物理，魔法，速度の最低ラインをそれぞれPin，Min，Sinとして，次の制約が増えます．

$Pin \leq P$
$Min \leq M$
$Sin \leq S$

以上の制約のもとで，今回は以下の目的関数を最大化します．

$f(a1,..,cn) = D = \sum_{i=1}^l ai * dai + \sum_{i=1}^m bi * dbi + \sum_{i=1}^n ci * dci$

あとはソルバーが勝手に解いてくれる

長くなりましたがやっと本題，以上で定式化した整数計画問題をPython-MIPで解いてみましょう．以下のコードで解くことができます．説明はコメントで書いています．

solver.py

from mip import Model, maximize, BINARY, xsum, OptimizationStatus

# ユーザが自分で決める物理，魔法，速度の下限値
pin = 25
min = 40
sin = 60

# ユーザのプレーヤーのレベル
lebel = 35

# 全装備のデータ
l = 10
m = 10
n = 10
# 物理
pa = [3, 4, 2, 5, 6, 3, 5, 6, 3, 2]
pb = [5, 7, 1, 10, 3, 4, 5, 6, 2, 2]
pc = [12, 45, 54, 43, 43, 5, 23, 34, 42, 43]
# 魔法
ma = [4, 6, 3, 8, 9, 1, 4, 3, 6, 7]
mb = [2, 4, 10, 4, 4, 4, 4, 1, 7, 8]
mc = [100, 9, 12, 43, 32, 5, 33, 24, 12, 93]
# 速度
sa = [12, 34, 32, 43, 54, 23, 32, 12, 18, 73]
sb = [32, 32, 12, 56, 64, 33, 43, 54, 64, 12]
sc = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
# 防御力
da = [43, 43, 23, 75, 52, 12, 32, 54, 32, 21]
db = [32, 32, 32, 54, 65, 54, 42, 75, 32, 42]
dc = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
# ランク
ra = [12, 13, 11, 19, 11, 9, 6, 10, 6, 12]
rb = [12, 18, 18, 15, 20, 9, 10, 12, 14, 18]
rc = [20, 12, 14, 17, 19, 18, 14, 13, 14, 12]

# まずモデルを定義していきます
solver = Model("最強装備シミュレータ(仮)")

# 変数を定義します．今回はa, b, cです．
# solverのadd_varメソッドを呼び出すと引数に指定した名前の変数がモデルの中に作成されます．
# その変数へのアドレスが返ってくるので，それをa，b，cという配列に保存しておきます．
# 変数は0か1しかとらないのでBINARYという種類の変数として定義しています．
# INTEGERとして定義して，後で制約で0以上1以下と追加しても問題としては同じです．
# CONTINUOUSとして定義すると実数値を取ることができ，混合整数計画問題が解けます
a = [solver.add_var('a{}'.format(i), var_type=BINARY) for i in range(l)]
b = [solver.add_var('b{}'.format(i), var_type=BINARY) for i in range(m)]
c = [solver.add_var('c{}'.format(i), var_type=BINARY) for i in range(n)]

# モデルに制約を追加していきます．add_constrメソッドを呼び出します．
# 総和はfor文で計算しても良いのですが，
# Python-MIPが提供するxsum関数で計算した方が高速なのでそちらを利用します．

# 頭，体，武器は一つずつしか装備できません
solver.add_constr(xsum(a[i] for i in range(l)) <= 1)
solver.add_constr(xsum(b[i] for i in range(m)) <= 1)
solver.add_constr(xsum(c[i] for i in range(n)) <= 1)

# 物理，魔法，速度，防御力，ランクの合計を表す変数を用意します
P = solver.add_var('P')
solver.add_constr(xsum(a[i] * pa[i] for i in range(l)) + xsum(b[i] * pb[i] for i in range(m)) + xsum(c[i] * pc[i] for i in range(n)) == P)
M = solver.add_var('M')
solver.add_constr(xsum(a[i] * ma[i] for i in range(l)) + xsum(b[i] * mb[i] for i in range(m)) + xsum(c[i] * mc[i] for i in range(n)) == M)
S = solver.add_var('S')
solver.add_constr(xsum(a[i] * sa[i] for i in range(l)) + xsum(b[i] * sb[i] for i in range(m)) + xsum(c[i] * sc[i] for i in range(n)) == S)
D = solver.add_var('D')
solver.add_constr(xsum(a[i] * da[i] for i in range(l)) + xsum(b[i] * db[i] for i in range(m)) + xsum(c[i] * dc[i] for i in range(n)) == D)
R = solver.add_var('R')
solver.add_constr(xsum(a[i] * ra[i] for i in range(l)) + xsum(b[i] * rb[i] for i in range(m)) + xsum(c[i] * rc[i] for i in range(n)) == R)

# ランクの総和はレベル以下です
solver.add_constr(R <= lebel)

# 物理，魔法，速度の下限値を設定します
solver.add_constr(pin <= P)
solver.add_constr(min <= M)
solver.add_constr(sin <= S)

# 最後に目的関数を設定します
# 今回はDの最大化なので以下のように書きます
solver.objective = maximize(D)

# 定式化したモデルを解くにはこのメソッドを呼び出します．
# 返り値は結果のステータスです
# 計算の過程がログとして出力されます
status = solver.optimize()

# 最適解がもとまったなら何番目の装備を使用したのかを出力しましょう
if status == OptimizationStatus.OPTIMAL:
    head = 0
    for i in range(l):
        if a[i].x >= 0.99:
            head = i
            print("頭装備は{}番目のものを使用".format(head + 1))
            break
    body = 0
    for i in range(m):
        if b[i].x >= 0.99:
            body = i
            print("体装備は{}番目のものを使用".format(body + 1))
            break
    hand = 0
    for i in range(m):
        if c[i].x >= 0.99:
            hand = i
            print("武器は{}番目のものを使用".format(hand + 1))
            break
    print("物理:{}，魔法:{}，速度:{}，防御力:{}，ランク:{}".format(P.x, M.x, S.x, D.x, R.x))
else:
    print("そのような装備の組み合わせはありません")

これを実行してみると...

$ python3 solver.py
頭装備は8番目のものを使用
体装備は8番目のものを使用
武器は10番目のものを使用
物理:55.0，魔法:97.0，速度:67.0，防御力:132.0，ランク:34.0

いけてそうですね！
これであなたもシミュレータ作者！

さいごに

簡単に解けた割にどんな問題でも解けそうな気がしませんか？これが整数計画問題として定式化して解く魅力です．

本当はヒューリスティックや整数計画問題と分枝限定法の関係や混合整数計画問題，整数計画問題，0-1整数計画問題や線形計画法，最大流で解ける差分制約式系の話をしたかったのですが長くなって疲れたのでやめました．もし興味があればググってみてください．長かったですがお付き合いくださりありがとうございました．

今回書いたコードは一応ここにあげておきます．
もっとPython-MIPで解くプログラム例が見たいという方は公式ドキュメントや電タクがブロック積み替え問題を解いているコード等を参考にしてください