こんにちは！しゅんです！
今回の記事は2部グラフの最小重み完全マッチングついて解説していきます！
マッチングとはグラフ理論という数学の分野で登場する用語で、現実世界の様々な分野で応用できます。
マッチングには色々な種類がありますが、その中でも今回は完全グラフじゃない部グラフの最小重み完全マッチングを求める問題をpythonで解いてみたいと思います！
2部グラフの最小重み完全マッチングについてはこちらの記事で解説しています！
完全2部グラフの最小重み完全マッチング問題の定式化はこちらの記事で解説しています！
左右の頂点数が異なる完全2部グラフの最小重みマッチング問題の定式化はこちらの記事で解説しています！
それではやっていきましょう！
普段は統計検定2級の記事を書いてたりします。
ぜひ他の記事も読んでみてください！
このブログの簡単な紹介はこちらに書いてあります。
興味があったら見てみてください。
このブログでは経営工学を勉強している現役理系大学生が、経営工学に関することを色々話していきます！
ぼくが経営工学を勉強している中で感じたことや、興味深かったことを皆さんと共有出来たら良いなと思っています。
そもそも経営工学とは何なのでしょうか。Wikipediaによると
経営工学（けいえいこうがく、英: engineering management）は、人・材料・装置・情報・エネルギーを総合したシステムの設計・改善・確立に関する活動である。そのシステムから得られる結果を明示し、予測し、評価するために、工学的な分析・設計の原理・方法とともに、数学、物理および社会科学の専門知識と経験を利用する。
引用元 : 経営工学 – Wikipedia
長々と書いてありますが、要は経営、経済の課題を理系的な観点から解決する学問です。
2パターンの定式化を考える
今回は2パターンの定式化を考えたいと思います。1つ目は「辺集合を中心に考える方法」で、2つ目は「辺の重みをめっちゃ大きくする方法」です。
この記事では「辺集合を中心に考える方法」について解説したいと思います！
「辺の重みをめっちゃ大きくする方法」についてはこちらの記事で解説しています！
現在制作中
今回解く問題
今回は以下の2部グラフの最小重み完全マッチングを求めたいと思います。
この問題を整数計画問題として定式化すると以下のようになります。今回はこの整数計画問題をpythonで実行してみたいと思います。
\(\min \;\;\; \sum\limits_{(i,j) \in E} w_{ij}x_{ij}\)
\(\;\text{s.t.}\;\;\;\sum\limits_{(i,j) \in E} x_{ij}=1 \;\;\;\; (i = 0,1,2)\)
\(\;\;\;\;\;\;\;\;\;\sum\limits_{(i,j) \in E} x_{ij} = 1 \;\;\;\; (j = 0,1,2)\)
\(\;\;\;\;\;\;\;\;\; x_{ij} \in \{0,1\} \;\;(i=0,1,2 \; j=0,1,2)\)
\(w_{ij}\) : 辺\((i,j)\)の重み
\(E\) : グラフの辺集合
この問題の定式化については以下の記事で詳しく解説しています！
コードの全貌と実行結果
まず最初にコードの全貌と実行結果を示してから、コードと結果について1つずつ解説していきたいと思います。コードの全貌は以下のようになります。
## 事前準備
!pip install pulp # pulpをインストール
from pulp import * # pulpをインポート
## パラメータの設定
v1 = [0,1,2] # 左側の頂点を設定
v2 = [0,1,2] # 右側の頂点を設定
E = [(0,0),(0,1),(1,0),(1,2),(2,0),(2,1),(2,2)] # 辺を設定
# 各辺の重みを設定
w = {(0,0):3,
(0,1):2,
(1,0):1,
(1,2):4,
(2,0):2,
(2,1):1,
(2,2):3}
## 整数計画問題の設定
prob = pulp.LpProblem(sense = LpMinimize) # 問題をMinimizeに設定
x = LpVariable.dict("x_%s_%s", E, cat="Binary") # 変数の設定
prob += lpSum(x[(i,j)] * w[(i,j)] for (i,j) in E) # 目的関数の設定
for i in v1:
prob += lpSum(x[(i,j)] for j in v2 if (i,j) in E) == 1 # 左側の頂点に関する制約条件を設定
for j in v2:
prob += lpSum(x[(i,j)] for i in v1 if (i,j) in E) == 1 # 右側の頂点に関する制約条件を設定
result = prob.solve() # 問題を解く
## 結果の表示
print(LpStatus[result]) # 得られた結果の状態を表示（Optimalなら最適解が得られている）
for (i,j) in E:
print(f"x{(i,j)} =",x[(i,j)].value()) # 最適解（各変数の値）を表示
print("最適値 :", prob.objective.value()) # 最適値を表示このコードを実行すると以下のような結果が得られました。
それでは1つずつ解説していきます。
コードを1つずつ解説する
事前準備
## 事前準備
!pip install pulp # pulpをインストール
from pulp import * # pulpをインポート問題を解くための事前準備としてpulpをインポートします。pulpはpythonで数理最適化の問題を解くために使うライブラリです。pulpをインストールしてからインポートしています。
パラメータの設定
## パラメータの設定
v1 = [0,1,2] # 左側の頂点を設定
v2 = [0,1,2] # 右側の頂点を設定
E = [(0,0),(0,1),(1,0),(1,2),(2,0),(2,1),(2,2)] # 辺を設定
# 各辺の重みを設定
w = {(0,0):3,
(0,1):2,
(1,0):1,
(1,2):4,
(2,0):2,
(2,1):1,
(2,2):3}次に各パラメータを設定していきます。今回は完全グラフじゃない2部グラフを考えているので頂点、辺、各辺の重みを設定する必要があります。
v1が左側の頂点、v2が右側の頂点を表しています。Eが辺集合を表しています。wは各辺の重さを辞書形式で表しています。例えばw[(0,1)]は2ですが、これは辺\((0,1)\)の重みが2であることを表しています。
整数計画問題の設定
## 整数計画問題の設定
prob = pulp.LpProblem(sense = LpMinimize) # 問題をMinimizeに設定
x = LpVariable.dict("x_%s_%s", E, cat="Binary") # 変数の設定
prob += lpSum(x[(i,j)] * w[(i,j)] for (i,j) in E) # 目的関数の設定
for i in v1:
prob += lpSum(x[(i,j)] for j in v2 if (i,j) in E) == 1 # 左側の頂点に関する制約条件を設定
for j in v2:
prob += lpSum(x[(i,j)] for i in v1 if (i,j) in E) == 1 # 右側の頂点に関する制約条件を設定
result = prob.solve() # 問題を解く
それではいよいよ整数計画問題をpythonで記述していきたいと思います。少し長いのでさらにいくつかのパートに分けて説明したいと思います。
問題をMinimizeに設定
prob = pulp.LpProblem(sense = LpMinimize) # 問題をMinimizeに設定この行では問題の種類を設定しています。pulp.LpProblemで数理最適化の問題を作ることができます。かっこの中のsense=の所で問題の種類を設定できます。今回は最小値を求めたいのでLpMinimizeと設定します。
最大値を求めたい場合はLpMaximizeとします。
0-1変数の設定
x = LpVariable.dict("x_%s_%s", E, cat="Binary") # 変数の設定
次に0-1変数\(x_{ij}\in\{0,1\}\)を設定しています。変数はLpVariableで設定できます。変数をたくさん作りたい場合は.dictを使うと便利です。かっこの中は左から順番に（”変数の名前” 添字の集合, cat=”変数の種類”）を表しています。
上の例で言うと名前はx、添字の集合は辺集合E、変数はBinary（0-1変数）と設定しています。つまりこのコードだけで
\(x_{00},x_{01},x_{10},x_{12},x_{20},x_{21},x_{22}\)
を作成することができます。”変数の名前”の所で%sが登場していますが、こうすることで辺集合\(E\)中の要素\((i,j)\)を\(x\)の添え字にすることができます。
Binary : 0-1変数
Integer : 整数変数
Continuous : 連続変数
目的関数の設定
prob += lpSum(x[(i,j)] * w[(i,j)] for (i,j) in E) # 目的関数の設定これは本筋ではありませんが、目的関数と制約条件はprob+=を最初に付けます。probは問題をMinimizeに設定したときにprob=としたのでそれに合わせます。もしMinimizeの設定の所で
problem = pulp.LpProblem(sense = LpMinimize) # 問題をMinimizeに設定と書いたら目的関数と制約条件の設定の所ではproblem+=と書く必要があります。
肝心の目的関数の設定についてですが、\(\sum\limits w_{ij}x_{ij}\)を表すにはlpSumが便利です。イメージ的にはlpSumが\(\sum\)の役割をしてくれると思ってください。
for (i,j) in Eと書くことで\(\sum\limits_{(i,j) \in E}\)を表現することができます。
制約条件の設定
for i in v1:
prob += lpSum(x[(i,j)] for j in v2 if (i,j) in E) == 1 # 左側の頂点に関する制約条件を設定
for j in v2:
prob += lpSum(x[(i,j)] for i in v1 if (i,j) in E) == 1 # 右側の頂点に関する制約条件を設定
次に制約条件を設定していきます。まず左側の頂点に関する制約
\(\sum\limits_{(i,j) \in E}x_{ij}=1 \;\; (i = 0,1,2)\)
を設定します。先ほどと同じように\(\sum\)はlpSumで表現します。かっこの中を(x[(i,j)] for j in v2 if (i,j) in E)とすることで
\(\sum\limits_{(i,j) \in E}x_{ij}\)
を表現することができます。（v2=[0,1,2]）
「for j in v2」だけだと\((i,j) \notin E\)の辺\((i,j)\)に対しても\(x_{ij}\)を足せと言ってしまいます。ところが\((i,j) \notin E\)の辺\((i,j)\)に関しては変数\(x_{ij}\)を生成していないのでエラーが出てしまいます。
「if (i,j) in E」を加えることで\((i,j) \in E\)を満たす辺\((i,j)\)だけを考えることができるため、うまく動きます。
これでが1になれば良いので==1としています。
= 1にするとエラーが出るので気を付けてください。
右側の頂点に関しても同じように設定します。
問題を解く
result = prob.solve() # 問題を解く
prob.solve()で問題を解くことができます。結果をresultに格納しておきます。
結果の表示
## 結果の表示
print(LpStatus[result]) # 得られた結果の状態を表示（Optimalなら最適解が得られている）
for (i,j) in E:
print(f"x{(i,j)} =",x[(i,j)].value()) # 最適解（各変数の値）を表示
print("最適値 :", prob.objective.value()) # 最適値を表示
最後に結果を表示するコードについて説明したいと思います。まずprint(LpStatus[result])で得られた結果の状態を表示します。もしOptimalと返されたらちゃんと最適値が求まっています。
次に最適解、すなわち各変数\(x_{ij}\)の値を表示するように設定します。解の値は.value()で得ることができます。辺集合\(E\)中の任意の辺\((i,j)\)に対して\(x_{ij}\)の値を表示したいので、for (i,j) in Eでiとjを回しています。
最後に最適値を表示するように設定します。最適値はobject.value()で得ることができます。
結果を1つずつ解釈する
上のコードを実行したら以下の結果が出力されたので、この結果を1つずつ解釈していきましょう。
1行目～5行目は無視して大丈夫です。（pulpがちゃんとインストールされたことを表しています。）
6行目は結果の状態を表しています。Optimalと表示されているのでちゃんと最適解が得られていることが分かりますね。
7行目～13行目では最適解（各変数の値）を表しています。これを見るとx(0,1)、x(1,0)、x(2,2)の3つだけが1でそれ以外は0になっていますね。
これはつまり辺\((0,1)\)、辺\((1,0)\)、辺\((2,2)\)を選ぶのが最小重み完全マッチングだと言うことを表しています。
13行目は最適値を表しています。これを見ると今回の問題の最適値は6であることが分かりますね。
おわりに
いかがでしたか。
今回の記事では完全グラフじゃない2部グラフの最小重み完全マッチングについて解説していきました。
今後もこのような数理最適化に関する記事を書いていきます！
最後までこの記事を読んでくれてありがとうございました。
