整数計画法で研究室配属

目的や経緯

授業でいくつかの研究室に学生を配属することになった．システマチックな方法で，できるだけ不満の少ない配属方法を求めたい．そこで整数計画法をつかうことにする．

Excel addinのソルバを使えば簡単と思ってたのだが，実行してみると「変数セルが多すぎます。」とのメッセージが．変数の数，最大200までらしいので，今回の問題（学生数100人，研究室数10）には使えない．
CPLEXとか独立したソルバを使おうかと考えたが，どうせ入力ファイルもプログラムで自動作成しないといけないので，今回は，Pythonのプログラムの中からソルバを動かす方針でいく．
具体的には，PuLPパッケージを利用．
学生は各研究室の希望順を回答済み（1が最高）．

import pulp

data = (
    ("名前", "唯野研", "蟻巣川研", "獅子成研"),
    ("一堂零", 1, 2, 3),
    ("冷越豪", 1, 2, 3),
    ("出瀬潔", 1, 3, 2),
    ("大間仁", 3, 1, 2),
    ("物星大", 3, 1, 2)
)

配属方法の方針

配属は以下のように2段階で決定．

方針1：最悪ケースの最善化

配属された研究室の希望順に関して，学生の中でもっとも下の順位だった学生について考える．公平のため，このような最悪の学生の配属研究室に対する希望順が，できるだけ上になることを最優先する．

したがって，100人の学生のうち99人が希望順1位の研究室になったとしても，残り一人が10位の研究室になるのであれば，100人全員が9位の研究室に配属される方を優先する．

方針2：全体の割り当て

最初の段階で求めた最良の最悪順位を守ったうえで，全員ができるだけ上の希望の研究室になるようにしたい．ただし，学生A, Bが希望順でそれぞれ1位, 4位の研究室に配属されるよりも，2位と3位の研究室に配属される方が良いものと考える．つまり，下の順位を改善することを，上の順位を改善することより優先する．

定式化

定数
$m$: 学生数, $n$: 研究室数, $a_{i,j}$: 学生$i (0\leq i \leq m-1)$の研究室$j (0\leq j \leq n-1)$に対する希望順位($1\leq a_{i,j} \leq n$), $size_{min}$: 研究室への配属人数の下限，$size_{max}$: 研究室への配属人数の上限

変数
0-1変数$x_{i,j}$ $(1 \leq i \leq m, 1 \leq j \leq n)$で，学生$i$が研究室$j$に割り当てられたことを表す．1なら割り当て．
変数$worst$ $(1 \leq j \leq n)$で最悪ケースでの配属研究室の希望順位を表す．

問題1：最悪ケースの最善化

目的関数: $worst$ →最小化
制約
$x_{i,1} + x_{i,2} + ... + x_{i, n} = 1$ 学生$i$は必ずどこかに配属になる
$x_{1,j} + x_{2,j} + ... + x_{m, j} \geq size_{min}$, $x_{1,j} + x_{2,j} + ... + x_{m, j} \leq size_{max}$ 研究室$j$の配属人数は$size_{min}$以上，$size_{max}$以下
$a_{i,1} x_{i,1} + a_{i,2} x_{i,2} + ... + a_{i,n} x_{i, n} \leq worst$ 学生$i$の配属研究室の希望順位は$worst$かそれ以上

問題2：全体の割り当て

問題1を解いて得られた最悪ケースの希望順位，つまり，最適解での$worst$の値を，$bound$とする．

目的関数: $\sum_i a_{i,1}^\alpha x_{i,1} + a_{i,2}^\alpha x_{i,2} + ... + a_{i,n}^\alpha x_{i, n}$ →最小化

$\alpha$は1以上の実数．$\alpha=1$なら，2位から1位への改善と，3位から2位への改善を，同じ重要性とみなすことになる．$\alpha>1$なら下位の改善を優先．たとえば，$\alpha=2$なら，2位から1位への改善は目的関数の値を$2^2 - 1^2 = 3$下げるが，3位から2位への改善は$3^2 - 2^2 = 5$下げる．

制約
$x_{i,1} + x_{i,2} + ... + x_{i, n} = 1$ 学生$i$は必ずどこかに配属になる
$x_{1,j} + x_{2,j} + ... + x_{m, j} \geq size_{min}$, $x_{1,j} + x_{2,j} + ... + x_{m, j} \leq size_{max}$ 研究室$j$の配属人数は$size_{min}$以上，$size_{max}$以下
$a_{i,1} x_{i,1} + a_{i,2} x_{i,2} + ... + a_{i,n} x_{i, n} \leq bound$ 学生$i$の配属研究室の希望順位は$bound$かそれ以上

コードと実行例

研究室には同じ人数を配属することにする．余りが出る場合は，配属人数の差を1以下にする．以下は上の続き．

NUM_MEMBERS = len(data)-1  # 学生の数
NUM_GROUPS = len(data[0])-1  # グループの数
MIN_GROUP_SIZE = NUM_MEMBERS//NUM_GROUPS  # グループの人数の最小値
# グループの人数の最大値 (最大値 - 最小値 <= 1)
if NUM_MEMBERS % NUM_GROUPS == 0:
    MAX_GROUP_SIZE = MIN_GROUP_SIZE
else:
    MAX_GROUP_SIZE = MIN_GROUP_SIZE + 1

# 問題1
problem1 = pulp.LpProblem('FindWorstCase', pulp.LpMinimize)

# 0-1変数x[i,j]を宣言 i:学生, j:グループ
x = {}
for i in range(0, NUM_MEMBERS):
    for j in range(0, NUM_GROUPS):
        x[i, j] = pulp.LpVariable("x({:},{:})".format(i, j), 0, 1, 'Binary')

# 整数変数worstを宣言　値は1以上，グループの数以下
worst = pulp.LpVariable('worst', 1, NUM_GROUPS, 'Integer')

# 制約条件A: 各学生は正確に一つのグループに割当
for i in range(0, NUM_MEMBERS):
    problem1 += sum(x[i, j] for j in range(0, NUM_GROUPS)) == 1

# 制約条件B: 各グループの人数は，最低値以上，最大値以下
for j in range(0, NUM_GROUPS):
    problem1 += sum(x[i, j] for i in range(0, NUM_MEMBERS)) >= MIN_GROUP_SIZE
    problem1 += sum(x[i, j] for i in range(0, NUM_MEMBERS)) <= MAX_GROUP_SIZE

# 制約条件C: 各学生の配属される研究室は，希望順でworstと同じかそれより若い
for i in range(0, NUM_MEMBERS):
    problem1 += sum(data[i+1][j+1] * x[i, j]
                    for j in range(0, NUM_GROUPS)) <= worst

# 目的関数：worstつまり希望順でのワーストケース
problem1 += worst

# 問題を解く
problem1.solve()
print("Worst:", worst.value())

# 問題2
bound = worst.value()  # 求めたワーストケースの場合を，希望順のboundとする
problem2 = pulp.LpProblem('AsssignAll', pulp.LpMinimize)

# 制約条件A: 各学生は正確に一つのグループに割当
for i in range(0, NUM_MEMBERS):
    problem2 += sum(x[i, j] for j in range(0, NUM_GROUPS)) == 1

# 制約条件B: 各グループの人数は，最低値以上，最大値以下
for j in range(0, NUM_GROUPS):
    problem2 += sum(x[i, j] for i in range(0, NUM_MEMBERS)) >= MIN_GROUP_SIZE
    problem2 += sum(x[i, j] for i in range(0, NUM_MEMBERS)) <= MAX_GROUP_SIZE

# 制約条件: 各学生の配属される研究室は，希望順でboundと同じかそれより上位
for i in range(0, NUM_MEMBERS):
    problem2 += sum(data[i+1][j+1] * x[i, j]
                    for j in range(0, NUM_GROUPS)) <= bound

# 目的関数
WEIGHT_EXPONENT = 1.2  # 重みづけのためのα
problem2 += sum((data[i+1][j+1] ** WEIGHT_EXPONENT) * x[i, j]
                for i in range(0, NUM_MEMBERS) for j in range(0, NUM_GROUPS))

# 問題を解く
problem2.solve()
print("名前", "研究室", "順位")
for i in range(0, NUM_MEMBERS):
    for j in range(0, NUM_GROUPS):
        if x[i, j].value() > 0:
            print(data[i+1][0], data[0][j+1], data[i+1][j+1])

実行結果

Worst: 2.0
名前 研究室 順位
一堂零 唯野研 1
冷越豪 唯野研 1
出瀬潔 獅子成研 2
大間仁 蟻巣川研 1
物星大 蟻巣川研 1

注意

• たいていの場合，$bound = n$として，問題2だけ解いても（目的関数と$worst$の値に関して）同じ結果が得られる．
• ortoolsを使った方が，幸せになれるかも．