割り当て問題の整数計画法による応用例(Python-MIPによるペアマッチング)

はじめに

以前にPython-MIPを使用してTSPを解きましたが、 今回は割り当て問題を整数計画問題として定式化を行い、解いてみようと思います。 また実際に具体的な実世界に近い問題を設定し解いてみます。実際に行ったことは以下になります。

• 割り当て問題を解く
• ペアマッチング例
• 一般割り当て問題を解く
• 複数ペアマッチング例

割り当て問題

まず二部グラフについての割り当て問題を考えます。 二部グラフはhttps://en.wikipedia.org/wiki/Bipartite_graph のように頂点の集合を、辺で直接つながっていない二つの頂点の集合に分けることのできるグラフです。 この二つの集合同士の割り当てはマッチング問題に相当します。重み付き二部マッチングは辺にコストが設定されている二部グラフに対して全ての頂点のペアを作成します。 この時、辺のコストの和が最小(or 最大)になるように辺を選択することで最適なマッチングを導出します。

例としては、クラスで男女同数の場合に男女間で相性やテストの点数差などでコストを設定し、全員のペアを作りながらコストを最小(できるだけ仲の良い男女をペアにする)する。 他の例としては、医者と患者の割り当てとして医者の得意分野と患者の症状でコストを設定し、治療の担当ペアを作成する問題などが挙げられます。 このようなペアの作成する手法はいくつかあり整数計画法でも解けますが、ここではもっと効率の良いハンガリアン法(Hungary)を使用します。他の手法としては最小費用流から求める方法があり辺の数が少なければこちらの方がかなり効率が良くなると思います。

このような問題を解く最も単純な例として、正規分布からランダムにコストを設定して二部グラフを作成し、コストが最小となるようにマッチングを構築します。 ハンガリアン法を使用するにはscipyのlinear_sum_assignmentがあります。頂点数 $N$ に対して計算量は $N^{3}$ でやってくれます。

    import numpy as np
from scipy.optimize import linear_sum_assignment
import matplotlib.pyplot as plt

n = 100
cost_mat = np.random.normal(50,10,(n,n))
mx, mn, mean = (cost_mat.max(), cost_mat.min(), cost_mat.mean())
print(mx, mn, mean)

row, col = linear_sum_assignment(cost_mat)

plt.hist([cost_mat.flatten(),cost_mat[row, col]], bins=20, log=True, label=["all", "match"])
plt.legend()
plt.show()
    

下の図は各コストに対するヒストグラムです。マッチに使用する辺のコストが小さいものを使用しており、row/colのunique数を確認するとうまくいっているようです。

ペアマッチング例

次に現実に近い例として男女のマッチングを考えます。 どのようなペアを作りたいかでコストの設定をしますが、今回は単純に「身長差」と設定します。身長は正規分布から作成しました。 身長差が小さければ小さいほどコストを小さく設定することで、トータルとして男女間の身長差の和が最小となるペアを作成します。 具体的にはコストを身長の二乗の差の平方根をとっています。ただし、女性の方が男性より身長が高い場合コストを2倍に設定するようにしました。

    n = 100
def abs_sp(x,y):
    if x - y >= 0:
        return np.sqrt(x**2 - y**2)
    else:
        return np.sqrt(y**2 - x**2)*2

male_heights = np.random.normal(170,8,n)
female_heights = np.random.normal(160,8,n)
mesh_male, mesh_female = np.meshgrid(male_heights, female_heights)
v_abs = np.vectorize(abs_sp)
cost_mat = v_abs(mesh_male, mesh_female)

mx, mn, mean = (cost_mat.max(), cost_mat.min(), cost_mat.mean())
print(mx, mn, mean)


row, col = linear_sum_assignment(cost_mat)

plt.hist([cost_mat.flatten(),cost_mat[row, col]], bins=20, log=True, label=["all", "match"])
plt.legend()
plt.show()
    

こちらもヒストグラムで確認すると、コストの小さい辺でマッチングを行っていることがわかります。 実際にマッチングペアを見てみます。

    male_heights_lis = male_heights.tolist()
female_heights_lis = female_heights.tolist()
for r, c in zip(row, col):
    print(f"male: No.{c:>2} {male_heights_lis[c]:.2f}cm  female: No.{r:>2} {female_heights_lis[r]:.2f} cm")
    

    male: No.40 175.70cm  female: No. 0 152.72 cm
male: No.67 160.55cm  female: No. 1 160.49 cm
male: No.19 156.72cm  female: No. 2 155.33 cm
male: No.78 164.65cm  female: No. 3 160.46 cm
male: No.25 163.86cm  female: No. 4 163.66 cm
male: No.31 164.74cm  female: No. 5 160.30 cm
male: No. 4 175.51cm  female: No. 6 153.27 cm
male: No. 0 169.13cm  female: No. 7 169.13 cm
...
    

ここでは身長差にしましたが、コスト関数の作成の仕方は多様です。 両者の希望値とステータスとのギャップなどをコストとして考慮すれば、多様な性質のペアを作成することができます。

一般割り当て問題

先ほどまでは1対1の割り当てでしたが、今度は 1対複数 または 複数対1 で解いてみます。 また二部グラフの頂点数は半分半分で同数でしたが、数が同数でなくとも解けるようにします。 これはNP-hardでハンガリアン法や最小費用流の手法は使えませんので、整数計画法(MIP)を使用します。

複数ペアマッチング例

今度は男女1対1ペアではなく、1対複数 で割り当てを行うことでペアの候補として一人につき複数のペアができるように構成します。 それぞれコストは先ほど同じく身長差とし、身長は正規分布から作成します。

それぞれの男女ペア(i,j)に対するコストを$c_{ij}$、男性iを女性jに割り当てるがどうかのバイナリ変数を$y_{ij}$、 女性iを男性jに割り当てるがどうかのバイナリ変数を$x_{ij}$として、一人につき$P$ペアを作ることにします。 すると以下のような制約式を設定できます。

1. $minimize \sum _{ij}c_{ij}x_{ij}+c_{ji}y_{ji}$
2. $\sum _{j}c_{ji}y_{ji}=P$
3. $\sum _{i}c_{ij}x_{ij}=P$
4. $x_{ij}=y_{ji}$

式2,3は各人についてP人を割り当てる式になります。また式4は男女ペアがどちらか一方通行にならないようにする制約になります。 実際のコードは以下になります。 男女ともに人数を変更でき(男m人、女n人)、1人につき任意ペアを構成することができますが、制約によっては解がない場合があります。 また、ライブラリはPython-MIPのインターフェースからCBC solverを使用しています。

    import mip
def solve_matching(costs):
    n = len(costs)
    m = len(costs[0])

    model = mip.Model()
    x = [[model.add_var(var_type=mip.BINARY) for j in range(m)]
        for i in range(n)]

    y = [[model.add_var(var_type=mip.BINARY) for j in range(n)]
        for i in range(m)]

    model.objective = mip.minimize(
                        mip.xsum(costs[i][j]*x[i][j] for i in range(m) for j in range(n))
                        +mip.xsum(costs[j][i]*y[i][j] for i in range(n) for j in range(m))
                    )

    for i in range(m):
        model += mip.xsum(y[i][j] for j in range(n)) == 3

    for i in range(n):
        model += mip.xsum(x[i][j] for j in range(m)) == 3

    for i in range(n):
        for j in range(m):
            model += y[j][i] == x[i][j]

    model.threads = -1
    #model.verbose = 0

    model.optimize(max_seconds=1000)

    print("Status:", model.status)
    print("objective", model.objective_value)

    return x, y
    

結果は以下になります。 $P=3$に設定しているので男性側からも、女性側からも丁度3ペア出来ているのが確認できると思います。

    Status: OptimizationStatus.OPTIMAL
objective 5189.693057692778

male  0: female:  3  7 17 
male  1: female:  7 13 15 
male  2: female:  4  6 11 
male  3: female:  1  4 12 
male  4: female:  1  6 11 
male  5: female:  8 16 18 
male  6: female:  2  8 10 
male  7: female:  0  5 14 
male  8: female: 16 18 19 
male  9: female:  3  7 13 
male 10: female:  4 12 14 
male 11: female:  2  8  9 
male 12: female:  0  5 15 
male 13: female:  2 10 17 
male 14: female: 13 14 15 
male 15: female:  0  5 12 
male 16: female:  1  6 11 
male 17: female:  9 16 19 
male 18: female:  9 18 19 
male 19: female:  3 10 17 
    

まとめ

状況に合わせて適切なコスト関数と定式化を行えば、それに対応した割り当てが可能になります。 今回例示した男女ペアマッチングだけではなく、仕事のリソースの割り当てや配送計画時の荷物配分などにも適用が可能です。 ただ問題の規模が大きくなったり複雑になれば、問題のモデルを単純化したり別の手法にも頼る必要が出てくるでしょう。