整数計画ソルバーが実行毎に異なる最適解を出力する謎を解け！

はじめに

この記事は数理最適化Advent Calender 2022の1日目の記事です．

複数の最適解を持つ最適化問題

最適化問題が複数の最適解を持つ場合には，整数計画ソルバーは最適解を１つ出力してくれます．ところが，Pythonから整数計画ソルバーを繰り返し実行するとその度に異なる最適解を出力することがあります．「複数の最適解を持つ問題なんて重箱の隅をつつくようなレアケースじゃないか？」と突っ込む方もいるかと思いますが，シフトスケジューリングや時間割作成など割当問題を雛形とする現実問題ではしばしば生じる状況です．
例えば，以下の簡単な割当問題を考えましょう．m人の作業員にn個の仕事を割り当てる問題を考えます．ここで，（１）仕事jはちょうど1人の作業員に割り当てる，（２）作業員iに2個以上の仕事を割り当てるものとします．ただし，仕事の数が少なく作業員に十分な数の仕事を割当てられない場合には，制約（２）の違反量を最小化します．（現実的な問題ではありませんが分かり易さのため簡単な例を挙げています）

$$\begin{array}{lll} \textnormal{最小化} & \displaystyle\sum_{i=1}^m s_i\\ \textnormal{条件} & \displaystyle\sum_{i=1}^m x_{ij} = 1, & j=1,\dots,n,\\ & \displaystyle\sum_{j=1}^n x_{ij} + s_i \ge 2, & i=1,\dots,m,\\ & x_{i,j} \in \{0,1\}, & i=1,\dots,m, j=1,\dots,n. \end{array}$$

Python-MIPを用いたプログラムと実行結果

Python-MIPを使ってこの整数計画問題を解いてみましょう．下記のプログラムでは，作業員5人(A,B,C,D,E)，仕事7個(a,b,c,d,e,f,g)の例を考えます．

# ライブラリのインポート
import itertools
from mip import *
# 入力データ
set_agent = {"A", "B", "C", "D", "E"}
set_job = {"a", "b", "c", "d", "e", "f", "g"}
# モデルの作成
mip_model = Model()
# 変数の設定
x = {}
for agent, job in itertools.product(set_agent, set_job):
    x[agent,job] = mip_model.add_var(var_type='B', name='x({},{})'.format(agent,job))
s = {}
for agent in set_agent:
    s[agent] = mip_model.add_var(var_type='C', name='s({})'.format(agent))
# 制約条件の設定
for job in set_job:
    mip_model.add_constr(
        xsum(x[agent,job] for agent in set_agent) == 1, name='AGENT({})'.format(job)
    )
for agent in set_agent:
    mip_model.add_constr(
        xsum(x[agent,job] for job in set_job) + s[agent] >= 2, name='JOB({})'.format(agent)
    )
# 目的関数の設定
mip_model.objective = minimize(xsum(s[agent] for agent in s))
# LPファイルの出力
mip_model.write("sample.lp")
# ソルバーの実行
mip_model.optimize()
# 計算結果の出力
for agent,job in sorted(x):
    if x[agent,job].x > 0.5:
        print(agent,job)

このプログラムを実行すると以下の結果が得られます．最適値3と最適解(A,a)(A,e)(C,f)(D,b)(D,g)(E,c)(E,d)が出力されます．

Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
A a
A e
C f
D b
D g
E c
E d

ここで再度プログラムを実行すると，さきほどと異なる最適解(A,g)(C,b)(C,f)(D,d)(D,e)(E,a)(E,c)が出力されます．このように実行するたびに異なる最適解が出力されて再現性が担保できなくなっています．これは大変に使い勝手が悪いです．

Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
A g
C b
C f
D d
D e
E a
E c

ここで，整数計画ソルバーが実際に解いている整数計画問題を確認しましょう．整数計画ソルバーをコマンドラインから実行する際には，下記のフォーマットのLPファイルを入力データとして与えます．上記のプログラムではLPファイルsample.lpに整数計画問題を書き出しています．1回目の実行時に作成されたLPファイルは以下の通りです．

sample.lp

Minimize
  s(E) + s(B) + s(A) + s(D) + s(C)
Subject To
 AGENT(e): x(E,e) + x(B,e) + x(A,e) + x(D,e) + x(C,e) = 1
 AGENT(c): x(E,c) + x(B,c) + x(A,c) + x(D,c) + x(C,c) = 1
 AGENT(d): x(E,d) + x(B,d) + x(A,d) + x(D,d) + x(C,d) = 1
 AGENT(a): x(E,a) + x(B,a) + x(A,a) + x(D,a) + x(C,a) = 1
 AGENT(g): x(E,g) + x(B,g) + x(A,g) + x(D,g) + x(C,g) = 1
 AGENT(f): x(E,f) + x(B,f) + x(A,f) + x(D,f) + x(C,f) = 1
 AGENT(b): x(E,b) + x(B,b) + x(A,b) + x(D,b) + x(C,b) = 1
 JOB(E): x(E,e) + x(E,c) + x(E,d) + x(E,a) + x(E,g) + x(E,f) + x(E,b)
   + s(E) >= 2
 JOB(B): x(B,e) + x(B,c) + x(B,d) + x(B,a) + x(B,g) + x(B,f) + x(B,b)
   + s(B) >= 2
 JOB(A): x(A,e) + x(A,c) + x(A,d) + x(A,a) + x(A,g) + x(A,f) + x(A,b)
   + s(A) >= 2
 JOB(D): x(D,e) + x(D,c) + x(D,d) + x(D,a) + x(D,g) + x(D,f) + x(D,b)
   + s(D) >= 2
 JOB(C): x(C,e) + x(C,c) + x(C,d) + x(C,a) + x(C,g) + x(C,f) + x(C,b)
   + s(C) >= 2
Bounds
Binaries
 x(E,e) x(E,c) x(E,d) x(E,a) x(E,g) x(E,f) x(E,b) x(B,e) x(B,c) x(B,d)
 x(B,a) x(B,g) x(B,f) x(B,b) x(A,e) x(A,c) x(A,d) x(A,a) x(A,g) x(A,f)
 x(A,b) x(D,e) x(D,c) x(D,d) x(D,a) x(D,g) x(D,f) x(D,b) x(C,e) x(C,c)
 x(C,d) x(C,a) x(C,g) x(C,f) x(C,b)
End

次に2回目の実行時に作成されたLPファイルを確認します．

sample.lp

Minimize
  s(D) + s(B) + s(C) + s(E) + s(A)
Subject To
 AGENT(b): x(D,b) + x(B,b) + x(C,b) + x(E,b) + x(A,b) = 1
 AGENT(e): x(D,e) + x(B,e) + x(C,e) + x(E,e) + x(A,e) = 1
 AGENT(d): x(D,d) + x(B,d) + x(C,d) + x(E,d) + x(A,d) = 1
 AGENT(f): x(D,f) + x(B,f) + x(C,f) + x(E,f) + x(A,f) = 1
 AGENT(a): x(D,a) + x(B,a) + x(C,a) + x(E,a) + x(A,a) = 1
 AGENT(g): x(D,g) + x(B,g) + x(C,g) + x(E,g) + x(A,g) = 1
 AGENT(c): x(D,c) + x(B,c) + x(C,c) + x(E,c) + x(A,c) = 1
 JOB(D): x(D,b) + x(D,e) + x(D,d) + x(D,f) + x(D,a) + x(D,g) + x(D,c)
   + s(D) >= 2
 JOB(B): x(B,b) + x(B,e) + x(B,d) + x(B,f) + x(B,a) + x(B,g) + x(B,c)
   + s(B) >= 2
 JOB(C): x(C,b) + x(C,e) + x(C,d) + x(C,f) + x(C,a) + x(C,g) + x(C,c)
   + s(C) >= 2
 JOB(E): x(E,b) + x(E,e) + x(E,d) + x(E,f) + x(E,a) + x(E,g) + x(E,c)
   + s(E) >= 2
 JOB(A): x(A,b) + x(A,e) + x(A,d) + x(A,f) + x(A,a) + x(A,g) + x(A,c)
   + s(A) >= 2
Bounds
Binaries
 x(D,b) x(D,e) x(D,d) x(D,f) x(D,a) x(D,g) x(D,c) x(B,b) x(B,e) x(B,d)
 x(B,f) x(B,a) x(B,g) x(B,c) x(C,b) x(C,e) x(C,d) x(C,f) x(C,a) x(C,g)
 x(C,c) x(E,b) x(E,e) x(E,d) x(E,f) x(E,a) x(E,g) x(E,c) x(A,b) x(A,e)
 x(A,d) x(A,f) x(A,a) x(A,g) x(A,c)
End

どうやら，変数や制約の出現順が入れ替わっていますね．これが異なる最適解が出力されてしまう原因です．多くの整数計画ソルバーは同じ問題でも入力データの変数や制約の出現順によりアルゴリズムの挙動が変化します．そのため，複数の最適解を持つ問題では実行するたびに異なる最適解が出力されてしまうことがあります．

実行毎に異なる最適解が出力される原因とその対策

では，なぜ実行するたびに変数や制約の出現順が入れ替わってしまうんでしょうか？じつは整数計画ソルバーではなくPythonに原因があることが分かりました．Pythonは辞書や集合などのデータ構造にハッシュを利用しています．Python3.3以降では内部で利用されるハッシュがランダム化されており，辞書や集合などのデータ構造を利用すると実行のたびに異なるハッシュテーブルが生成されます．その結果，プログラムを実行するたびに変数や制約を設定する際の処理順も入れ替わってしまったのです．

記事を公開後にご指摘をいただきました．Python3.7から辞書の順番は保持されるようになったそうです．自身で確認した範囲では大丈夫なようでした．

これを解決する方法は２つあります．
１つ目はプログラムを実行する前に環境変数PYTHONHASHSEEDの値を設定し，ハッシュのランダム化に使われている乱数の種を固定することです．プログラムを実行する環境によりますが，例えばMacのターミナルであれば以下のように設定します．

$ export PYTHONHASHSEED=0

ただし，プログラム実行前に環境変数PYTHONHASHSEEDの値をいちいち設定する必要があり面倒です．プログラムの先頭に下記のコードを追加すれば良いのでは？と思うかも知れませんが残念ながら上手く行きません．

import os
os.environ['PYTHONHASHSEED'] = '0'

2つ目は変数・制約・目的関数を設定する際に，for文などの繰り返し処理で辞書や集合を整列することです．このように，処理順を固定すれば変数や制約の出現順を固定できます．1つ目の方法に比べると修正する箇所は増えますがプログラムの修正だけで済みます．

# 変数の設定
x = {}
for agent, job in sorted(itertools.product(set_agent, set_job)):
    x[agent,job] = mip_model.add_var(var_type='B', name='x({},{})'.format(agent,job))
s = {}
for agent in sorted(set_agent):
    s[agent] = mip_model.add_var(var_type='C', name='s({})'.format(agent))
# 制約条件の設定
for job in sorted(set_job):
    mip_model.add_constr(
        xsum(x[agent,job] for agent in set_agent) == 1, name='AGENT({})'.format(job)
    )
for agent in sorted(set_agent):
    mip_model.add_constr(
        xsum(x[agent,job] for job in set_job) + s[agent] >= 2, name='JOB({})'.format(agent)
    )
# 目的関数の設定
mip_model.objective = minimize(xsum(s[agent] for agent in sorted(s)))

ちなみに，今回のプログラムでは作業員と仕事の集合を下記のようにリストで定義すれば繰返し実行しても同じ最適解を出力してくれそうですが，現実問題を整数計画ソルバーで解くのにリストしか使えないのは不便な話です．（確実な方法なのかも分かりません）

# 入力データ
set_agent = ["A", "B", "C", "D", "E"]
set_job = ["a", "b", "c", "d", "e", "f", "g"]

おわりに

本記事では，整数計画ソルバーが実行毎に異なる最適解を出力してしまう原因とその対策を紹介しました．私自身はずっと整数計画ソルバーの問題だと思っていたのですが，原因は整数計画ソルバーではなくPythonにあったというオチです．そのような次第で，整数計画ソルバーに限らず辞書や集合などPythonが内部でハッシュを利用するデータ構造では再現性が担保されない場合があることに注意してください．