●ハミルトン路
あるグラフにおいて、すべての頂点をちょうど一回ずつ通る経路を「ハミルトン路」といいます。また、すべての頂点を一回ずつ通って出発点に戻る経路を「ハミルトン閉路」といいます。任意のグラフがハミルトン路を持っているか否かを判定する方法は、現在では知られていません (NP 問題といわれています)。
今回は以下の経路図で A から I までのハミルトン路を求めてみましょう。
Ａ───Ｂ───Ｃ
│ │ │
│ │ │
│ │ │
Ｄ───Ｅ───Ｆ
│ │ │
│ │ │
│ │ │
Ｇ───Ｈ───Ｉ
図 : 経路図
●プログラムの作成
拙作のページ「最短経路問題」のプログラムは、スタートとゴール以外の頂点の制約条件で、入力の辺と出力の辺が等しくなることを指定しました。今回はすべての頂点を通るので、入力の辺と出力の辺が両方とも 1 になることを条件に指定すれば解けるような気がします。実際にプログラムを作って確かめてみましょう。
リスト : ハミルトン路を求める (No Good)
import pulp
#
# 隣接行列
# A B C D E F G H I
ws = [[ 0, 1, 0, 1, 0, 0, 0, 0, 0], # A
[ 1, 0, 1, 0, 1, 0, 0, 0, 0], # B
[ 0, 1, 0, 0, 0, 1, 0, 0, 0], # C
[ 1, 0, 0, 0, 1, 0, 1, 0, 0], # D
[ 0, 1, 0, 1, 0, 1, 0, 1, 0], # E
[ 0, 0, 1, 0, 1, 0, 0, 0, 1], # F
[ 0, 0, 0, 1, 0, 0, 0, 1, 0], # G
[ 0, 0, 0, 0, 1, 0, 1, 0, 1], # H
[ 0, 0, 0, 0, 0, 1, 0, 1, 0], # I
]
node = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I']
size = len(node) # 頂点の数
start = 0 # スタート (A)
goal = 8 # ゴール (I)
prob = pulp.LpProblem('hamilton')
# 変数の生成
xs = []
for i in range(size):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size)]
xs.append(xs1)
# 制約条件
prob += pulp.lpSum(xs[start]) == 1
prob += pulp.lpSum([xs[i][goal] for i in range(size)]) == 1
for i in range(size):
if i == start or i == goal: continue
prob += pulp.lpSum([xs[j][i] for j in range(size)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size)]) == 1
# 重みが 0 の辺は選択しない
for i in range(size):
for j in range(size):
if ws[i][j] == 0: prob += xs[i][j] == 0
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
if xs[i][j].value() == 1:
print('{} -> {}'.format(node[i], node[j]))
実行結果は次のようになりました。
Status Optimal A -> B B -> C C -> F D -> G E -> H F -> I G -> D H -> E
A -> B -> C -> F -> I という経路と、2 つの閉路 D <-> G, E <-> H ができるため、ハミルトン路にはなりません。
A -- B -- C
|
D E F
| | |
G H I
このように、頂点に出入りする辺の本数だけでは条件不足になるのです。
●閉路のチェック
閉路ができるのを防ぐため、今回は経路に沿って頂点に番号を順番に付けていくことにします。このとき、隣り合った頂点の差分は 1 になることに注意してください。閉路の場合、たとえば B - C - F - E で番号を 2, 3, 4, 5 と付けると、B と E は隣り合っているにもかかわらず差分が 3 になるため、この経路は除外されることになります。
部分巡回路を取り除くための制約条件を「部分巡回路除去制約」といいます。拙作のページ「制約論理プログラミング超入門: 経路の探索 (ハミルトン路)」ではこの方法でうまくいきましたが、PuLP (CBC) ではこれをどうやってプログラムしたらよいのか、皆目見当がつきません。いろいろと検索したところ、参考 URL 5, 6, 7 (pulp01.html#cite) に部分巡回路除去制約の説明がありました。有益な情報を公開されている作者の皆様さんに感謝いたします。
詳細は上記参考 URL をお読みいただくとして、ここでは制約条件を簡単に説明します。頂点 i, j に対応する変数を \(U_i \ U_j\) とします。頂点数が M 個、変数の値が 0 から M - 1 までとすると、\(U_i\) と \(U_j\) の間には次の関係が成り立ちます。
\(
U_i + 1 - bigM \times (1 - X_{ij}) \leq U_j
\)
これを「ポテンシャル制約」とか「Miller-Tucker-Zemlin (MTZ) 方式」と呼びます。\(X_{ij}\) が 1 の場合 (i から j につながる辺を選ぶ場合)、上式は \(U_i + 1 \leq U_j\) となります。これは変数の値が経路に沿って昇順に並ぶことを表します。
次に、\(X_{ij}\) が 0 の場合ですが、上式は \(U_i + 1 - bigM \leq U_j\) になります。変数の値は \(U_0\) を 0 とすると、それ以外の変数の値は \(1 \leq U_i \leq M - 1 \ (i \gt 0)\) になります。ここで、bigM の値を M - 1 とすると、左辺がとりうる値は 3 - M から 1 までになります。つまり、\(X_{ij}\) が 0 であれば上式は必ず成立します。
あとは \(U_i\) と \(U_j\) の間に上式の制約条件を設定すれば、経路に沿って頂点に番号を順番に付けていくことができます。ただし、参考 URL 6 によると MTZ 方式の制約は弱いことが知られているそうです。今回は簡単な例題なので、上式をそのまま使うことにします。
プログラムは次のようになります。
リスト : ハミルトン路を求める (Good)
import pulp
#
# 隣接行列
# A B C D E F G H I
ws = [[ 0, 1, 0, 1, 0, 0, 0, 0, 0], # A
[ 1, 0, 1, 0, 1, 0, 0, 0, 0], # B
[ 0, 1, 0, 0, 0, 1, 0, 0, 0], # C
[ 1, 0, 0, 0, 1, 0, 1, 0, 0], # D
[ 0, 1, 0, 1, 0, 1, 0, 1, 0], # E
[ 0, 0, 1, 0, 1, 0, 0, 0, 1], # F
[ 0, 0, 0, 1, 0, 0, 0, 1, 0], # G
[ 0, 0, 0, 0, 1, 0, 1, 0, 1], # H
[ 0, 0, 0, 0, 0, 1, 0, 1, 0], # I
]
node = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I']
size = len(node) # 頂点の数
start = 0 # スタート (A)
goal = 8 # ゴール (I)
prob = pulp.LpProblem('hamilton1')
# 変数の生成
xs = []
for i in range(size):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size)]
xs.append(xs1)
# 頂点に対応する変数
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size-1) for i in range(size)]
# 制約条件
prob += pulp.lpSum(xs[start]) == 1
prob += pulp.lpSum([xs[i][goal] for i in range(size)]) == 1
for i in range(size):
if i == start or i == goal: continue
prob += pulp.lpSum([xs[j][i] for j in range(size)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size)]) == 1
# 重み 0 の辺は選択しない
for i in range(size):
for j in range(size):
if ws[i][j] == 0: prob += xs[i][j] == 0
# 頂点に対応する変数の制約
for i in range(size):
for j in range(1, size):
if i == j: continue
prob += us[i] + 1 - (size - 1) * (1 - xs[i][j]) <= us[j]
prob += us[0] == 0
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
if xs[i][j].value() == 1:
print('{} -> {}'.format(node[i], node[j]))
print([u.value() for u in us])
頂点に対応する変数 (Integer) は配列 us に格納します。あとは、us の制約条件を設定するだけですが、j が 0 のときの制約は省くことに注意してください。MTZ 方式の左辺は最大で 1 になりますが、右辺が 0 (us[0] == 0) になると条件が不成立になるため、問題を解くことができなくなります。あとは、難しいところはないと思います。
それでは実行してみましょう。
Status Optimal A -> B B -> C C -> F D -> G E -> D F -> E G -> H H -> I [0.0, 1.0, 2.0, 5.0, 4.0, 3.0, 6.0, 7.0, 8.0]
A -- B -- C
|
D -- E -- F
|
G -- H -- I
今度はハミルトン路を求めることができました。us の値は経路 A - B - C - F - E - D - G - H - I に沿って、順番に番号が割り当てられていることがわかります。
●ハミルトン閉路
次は下図に示す経路図でハミルトン閉路を求めてみましょう。
Ａ───Ｂ───Ｃ │ │ ／│ │ │ ／ │ │ │／ │ Ｄ───Ｅ───Ｆ │ ／│ │ │ ／ │ │ │／ │ │ Ｇ───Ｈ───Ｉ 図 : 経路図
ハミルトン閉路のプログラムは、ハミルトン路のプログラムを改造すると簡単です。今回求める経路は閉路なので、すべての頂点に辺が 2 本つながることになります。これは、スタートとゴールの制約を削除すれば OK です。最後に、頂点 (M - 1) から最初の頂点 (0) につながる辺を選択して閉路は完成しますが、頂点 0 に戻る辺の制約条件は設定していないので、このままで大丈夫です。
プログラムは次のようになります。
リスト : ハミルトン閉路を求める
import pulp, time
#
# 隣接行列
# A B C D E F G H I
ws = [[0, 1, 0, 1, 0, 0, 0, 0, 0], # A
[1, 0, 1, 0, 1, 0, 0, 0, 0], # B
[0, 1, 0, 0, 1, 1, 0, 0, 0], # C
[1, 0, 0, 0, 1, 0, 1, 0, 0], # D
[0, 1, 1, 1, 0, 1, 1, 1, 0], # E
[0, 0, 1, 0, 1, 0, 0, 0, 1], # F
[0, 0, 0, 1, 1, 0, 0, 1, 0], # G
[0, 0, 0, 0, 1, 0, 1, 0, 1], # H
[0, 0, 0, 0, 0, 1, 0, 1, 0], # I
]
node = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I']
size = len(node) # 頂点の数
prob = pulp.LpProblem('hamilton2')
# 変数の生成
xs = []
for i in range(size):
xs1 = [pulp.LpVariable('v{}_{}'.format(i, j), cat='Binary') for j in range(size)]
xs.append(xs1)
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size-1) for i in range(size)]
# 制約条件
for i in range(size):
prob += pulp.lpSum([xs[j][i] for j in range(size)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size)]) == 1
# 重み 0 の辺は選択しない
for i in range(size):
for j in range(size):
if ws[i][j] == 0: prob += xs[i][j] == 0
# 頂点に対応する変数の制約
for i in range(size):
for j in range(1, size):
if i == j: continue
prob += us[i] + 1 - (size - 1) * (1 - xs[i][j]) <= us[j]
prob += us[0] == 0
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
if xs[i][j].value() == 1:
print('{} -> {}'.format(node[i], node[j]))
print([u.value() for u in us])
実行結果は次のようになりました。
Status Optimal A -> D B -> A C -> B D -> G E -> H F -> C G -> E H -> I I -> F [0.0, 8.0, 7.0, 1.0, 3.0, 6.0, 2.0, 4.0, 5.0]
A--B--C | | | | D E F | /| | |/ | | G H--I
PuLP (CBC) でもハミルトン閉路を求めることができました。
●騎士の巡歴
ナイト (騎士) はチェスの駒のひとつで将棋の桂馬の動きを前後左右にとることができます。次の図を見てください。
┌─┬─┬─┬─┬─┐ ┌─┬─┬─┬─┬─┐
│ │●│ │●│ │ │Ｓ│ │ │ │ │
├─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┤
│●│ │ │ │●│ │ │ │ │ │ │
├─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┤
│ │ │Ｋ│ │ │ │ │ │ │ │ │
├─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┤
│●│ │ │ │●│ │ │ │ │ │ │
├─┼─┼─┼─┼─┤ ├─┼─┼─┼─┼─┤
│ │●│ │●│ │ │ │ │ │ │Ｇ│
└─┴─┴─┴─┴─┘ └─┴─┴─┴─┴─┘
●：ナイト (K) が動ける位置 問題
問題 : 騎士の巡歴
「騎士の巡歴 (Knight Tour)」は、ナイトを動かして N 行 N 列の盤面のどのマスにもちょうど一回ずつ訪れるような経路を求める問題です。ちなみに、3 行 3 列、4 行 4 列の盤面には解がありませんが、5 行 5 列の盤面には解があります。今回は条件をひとつ追加して、スタート (S) からゴール (G) までの経路を求めることにします。この場合、N が偶数だと解はありません。これは簡単に証明できるので、興味のある方は考えてみてください。
●隣接行列の作成
┌─┬─┬─┬─┬─┐ │０│１│２│３│４│ ├─┼─┼─┼─┼─┤ │５│６│７│８│９│ ├─┼─┼─┼─┼─┤ │10│11│12│13│14│ ├─┼─┼─┼─┼─┤ │15│16│17│18│19│ ├─┼─┼─┼─┼─┤ │20│21│22│23│24│ 図 : マスの番号 └─┴─┴─┴─┴─┘
最初に、N 行 N 列盤の隣接行列 (N2 * N2 行列) を生成するプログラムを作りましょう。たとえば、5 行 5 列盤の場合、上図のようにマスに番号を付けることにします。
この場合、0 から移動できる場所は 7 と 11 に、12 から移動できる場所は 1, 3, 5, 9, 15, 19, 21, 23 になります。
プログラムは次のようになります。
リスト : 隣接行列の生成
def make_matrix(size):
size2 = size * size
ws = [[0] * size2 for _ in range(size2)]
for y in range(size):
for x in range(size):
k = y * size + x
for dx, dy in [(-1, -2), (1, -2), (-2, -1), (2, -1),
(-2, 1), (2, 1), (-1, 2), (1, 2)]:
x1 = x + dx
y1 = y + dy
if 0 <= x1 < size and 0 <= y1 < size:
ws[k][y1 * size + x1] = 1
return ws
関数 make_matrix の引数 size が盤面の大きさ (size * size) を表します。size2 = size * size とすると、隣接行列 ws の大きさは size2 * size2 になります。ws の要素は 0 で初期化します。あとは、盤面の座標 x, y と騎士の移動方向 dx, dy から移動先の座標 x1, y1 を求め、それが盤面の範囲内に収まっていれば、その位置に対応する ws の値を 1 にセットします。最後に ws を返します。
簡単な実行例を示しましょう。
>>> import knight >>> for x in knight.make_matrix(3): print(x) ... [0, 0, 0, 0, 0, 1, 0, 1, 0] [0, 0, 0, 0, 0, 0, 1, 0, 1] [0, 0, 0, 1, 0, 0, 0, 1, 0] [0, 0, 1, 0, 0, 0, 0, 0, 1] [0, 0, 0, 0, 0, 0, 0, 0, 0] [1, 0, 0, 0, 0, 0, 1, 0, 0] [0, 1, 0, 0, 0, 1, 0, 0, 0] [1, 0, 1, 0, 0, 0, 0, 0, 0] [0, 1, 0, 1, 0, 0, 0, 0, 0] >>> for x in knight.make_matrix(4): print(x) ... [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0] [0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0] [0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0] [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0] [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0] [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0] [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1] [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0] [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0] [1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1] [0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0] [0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0] [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0] [0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0] [0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0] [0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0]
●解法プログラム
隣接行列ができると、あとのプログラムは拙作のページ「ハミルトン路」で作成したプログラムとほとんど同じです。解法プログラムは次のようになります。
リスト ; 騎士の巡歴
def solver(size):
size2 = size * size
start = 0 # 右上隅
goal = size2 - 1 # 左下隅
ws = make_matrix(size) # 隣接行列
prob = pulp.LpProblem('knight-tour')
# 変数の生成
xs = []
for i in range(size2):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size2)]
xs.append(xs1)
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size2-1) for i in range(size2)]
# 制約条件
prob += pulp.lpSum(xs[start]) == 1
prob += pulp.lpSum([xs[i][goal] for i in range(size2)]) == 1
for i in range(size2):
if i == start or i == goal: continue
prob += pulp.lpSum([xs[j][i] for j in range(size2)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size2)]) == 1
# 重み 0 の辺は選択しない
for i in range(size2):
for j in range(size2):
if ws[i][j] == 0: prob += xs[i][j] == 0
# 頂点に対応する変数の制約の制約
for i in range(size2):
for j in range(1, size2):
if i == j: continue
prob += us[i] + 1 - (size2 - 1) * (1 - xs[i][j]) <= us[j]
for i in range(1, size2): 1 <= us[i]
prob += us[0] == 0 # スタート
prob += us[size2 - 1] == size2 - 1 # ゴール
prob += xs[0][size * 2 + 1] == 1 # 最初に選ぶ辺を限定
prob += us[size * 2 + 1] == 1
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
k = i * size + j
print('{:3d}'.format(int(us[k].value())), end=" ")
print("")
スタートから出る辺は二つありますが、それを一つに限定しています。それから、制約条件 \(1 \leq us[i] \ (i \ne 0)\) も追加しています。これがないと実行時間が少し遅くなるようです。
●実行結果
実行結果は次のようになりました。
>>> import knight1 >>> import time >>> s = time.time(); knight1.solver(5); time.time() - s Status Optimal 0 19 10 5 22 9 4 21 18 11 20 1 14 23 6 15 8 3 12 17 2 13 16 7 24 0.8754072189331055 >>> s = time.time(); knight1.solver(7); time.time() - s Status Optimal 0 27 2 5 14 25 10 3 16 7 26 11 40 13 28 1 4 15 6 9 24 17 36 29 8 41 12 39 30 33 46 37 20 23 42 35 18 31 44 47 38 21 32 45 34 19 22 43 48 1.5901174545288086 >>> s = time.time(); knight1.solver(9); time.time() - s Status Optimal 0 5 68 49 70 73 10 15 18 67 50 71 4 9 16 19 74 11 6 1 48 69 72 59 12 17 14 51 66 3 8 45 20 75 60 33 2 7 44 47 58 39 34 13 76 43 52 65 28 21 46 77 32 61 24 27 22 55 40 57 38 35 78 53 42 25 64 29 36 79 62 31 26 23 54 41 56 63 30 37 80 1.4380414485931396 >>> s = time.time(); knight1.solver(11); time.time() - s Status Optimal 0 99 90 5 28 93 54 7 50 45 48 89 4 101 94 91 6 29 44 47 8 51 100 1 98 27 102 55 92 53 112 49 46 97 88 3 32 95 26 103 30 43 52 9 2 71 96 87 56 31 42 111 10 113 116 73 86 33 68 41 110 25 104 115 12 23 70 67 72 57 34 105 40 11 24 117 114 85 74 69 80 109 58 35 118 39 22 13 66 77 108 83 106 81 62 17 20 119 38 75 84 79 64 59 16 19 36 61 14 21 78 65 76 107 82 63 60 15 18 37 120 2.028217077255249 >>> s = time.time(); knight1.solver(13); time.time() - s Status Optimal 0 139 10 131 6 141 8 151 78 157 76 161 80 121 132 3 140 9 130 5 156 145 150 79 158 75 138 1 122 11 4 7 142 149 152 77 160 81 162 123 120 133 2 129 38 61 144 155 146 153 74 159 134 137 124 127 12 143 148 39 62 69 84 163 82 119 126 135 98 37 128 13 60 147 154 63 70 73 136 97 118 125 100 35 40 105 68 85 72 83 164 117 108 99 36 41 106 101 14 59 18 67 64 71 112 25 96 107 102 15 34 17 104 57 86 165 66 109 116 111 26 95 42 103 58 19 166 65 48 87 24 113 28 31 92 33 16 43 54 89 56 167 50 29 110 115 22 27 94 91 20 45 52 49 88 47 114 23 30 93 32 21 44 53 90 55 46 51 168 24.058934688568115
- 実行環境 : Python 3.10.6, PuLP 2.7.0, Ubunts 22.04 LTS (WSL2, Windows 10), Intel Core i5-6200U 2.30GHz
N が大きくなると時間がかかるようになりますが、N = 13 でも 24 秒で解を求めることができました。
ご参考までに、Python/Tkinter による描画プログラムを示します。
リスト : Tkinter による経路の描画
import tkinter as tk
def draw(size, xs):
root = tk.Tk()
c0 = tk.Canvas(root, width = 40 * size + 20, height = 40 * size + 20)
c0.pack()
for y in range(size):
for x in range(size):
if (y + x) % 2 == 0:
color = 'white'
else:
color = 'black'
x1 = x * 40 + 10
y1 = y * 40 + 10
c0.create_rectangle(x1, y1, x1 + 40, y1 + 40, fill = color)
for i in range(size ** 2):
for j in range(size ** 2):
if xs[i][j].value() == 0: continue
y1 = (i // size) * 40 + 30
x1 = (i % size) * 40 + 30
y2 = (j // size) * 40 + 30
x2 = (j % size) * 40 + 30
c0.create_line(x1, y1, x2, y2, fill = 'pink', width = 3.0)
root.mainloop()
13 行 13 列盤の解の一例
●プログラムリスト
#
# knight1.py : 騎士の巡歴
#
# Copyright (C) 2019-2023 Makoto Hiroi
#
import pulp
import tkinter as tk
# 描画
def draw(size, xs):
root = tk.Tk()
c0 = tk.Canvas(root, width = 40 * size + 20, height = 40 * size + 20)
c0.pack()
for y in range(size):
for x in range(size):
if (y + x) % 2 == 0:
color = 'white'
else:
color = 'black'
x1 = x * 40 + 10
y1 = y * 40 + 10
c0.create_rectangle(x1, y1, x1 + 40, y1 + 40, fill = color)
for i in range(size ** 2):
for j in range(size ** 2):
if xs[i][j].value() == 0: continue
y1 = (i // size) * 40 + 30
x1 = (i % size) * 40 + 30
y2 = (j // size) * 40 + 30
x2 = (j % size) * 40 + 30
c0.create_line(x1, y1, x2, y2, fill = 'pink', width = 3.0)
root.mainloop()
# 隣接行列の生成
def make_matrix(size):
size2 = size * size
ws = [[0] * size2 for _ in range(size2)]
for y in range(size):
for x in range(size):
k = y * size + x
for dx, dy in [(-1, -2), (1, -2), (-2, -1), (2, -1),
(-2, 1), (2, 1), (-1, 2), (1, 2)]:
x1 = x + dx
y1 = y + dy
if 0 <= x1 < size and 0 <= y1 < size:
ws[k][y1 * size + x1] = 1
return ws
def solver(size):
size2 = size * size
start = 0 # 右上隅
goal = size2 - 1 # 左下隅
ws = make_matrix(size) # 隣接行列
prob = pulp.LpProblem('knight-tour')
# 変数の生成
xs = []
for i in range(size2):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size2)]
xs.append(xs1)
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size2-1) for i in range(size2)]
# 制約条件
prob += pulp.lpSum(xs[start]) == 1
prob += pulp.lpSum([xs[i][goal] for i in range(size2)]) == 1
for i in range(size2):
if i == start or i == goal: continue
prob += pulp.lpSum([xs[j][i] for j in range(size2)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size2)]) == 1
# 重み 0 の辺は選択しない
for i in range(size2):
for j in range(size2):
if ws[i][j] == 0: prob += xs[i][j] == 0
# 頂点に対応する変数の制約の制約
for i in range(size2):
for j in range(1, size2):
if i == j: continue
prob += us[i] + 1 - (size2 - 1) * (1 - xs[i][j]) <= us[j]
for i in range(1, size2): 1 <= us[i]
prob += us[0] == 0 # スタート
prob += us[size2 - 1] == size2 - 1 # ゴール
prob += xs[0][size * 2 + 1] == 1 # 最初に選ぶ辺を限定
prob += us[size * 2 + 1] == 1
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
k = i * size + j
print('{:3d}'.format(int(us[k].value())), end=" ")
print("")
# Tkinter による描画
# draw(size, xs)
●騎士の周遊
ところで、騎士の巡歴は「どのマスにもちょうど一回ずつ訪れたのち最初のマスに戻ってくること」を条件にする場合があります。これを「騎士の周遊」と呼びます。この場合、4 行 4 列盤や 5 行 5 列盤には解がありません。また、N 行 N 列の盤面で N が奇数の場合も、騎士は出発点に戻ることはできません。これも簡単に証明することができます。興味のある方は考えてみてください。
今回は「騎士の周遊」を解くプログラムを作りましょう。
●プログラムの作成
基本的には拙作のページ「ハミルトン閉路」と同じですが、スタート地点の制約条件を追加します。左上隅 (0 番目) をスタート (us[0] == 0) とすると、ここには 2 つの辺しかないので、閉路を生成するにはどちらの辺も選ぶことになります。そして、頂点の番号もどちらかが 1 で他方が size * size - 1 になります。実際には最初に選ぶ辺を限定して制約条件を追加します。
プログラムは次のようになります。
リスト : 騎士の周遊
def solver(size):
size2 = size * size
ws = make_matrix(size) # 隣接行列
prob = pulp.LpProblem('knight-circle')
# 変数の生成
xs = []
for i in range(size2):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size2)]
xs.append(xs1)
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size2-1) for i in range(size2)]
# 制約条件
for i in range(size2):
prob += pulp.lpSum([xs[j][i] for j in range(size2)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size2)]) == 1
# 重み 0 の辺は選択しない
for i in range(size2):
for j in range(size2):
if ws[i][j] == 0: prob += xs[i][j] == 0
for i in range(size2):
for j in range(1, size2):
if i == j: continue
prob += us[i] + 1 - (size2 - 1) * (1 - xs[i][j]) <= us[j]
prob += us[0] == 0 # スタート
for i in range(1, size2): 1 <= us[i]
prob += xs[0][size * 2 + 1] == 1 # 最初に選ぶ辺を限定
prob += us[size * 2 + 1] == 1
prob += xs[size + 2][0] == 1 # 最後に選ぶ辺
prob += us[size + 2] == size2 - 1 # 最後の頂点
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
k = i * size + j
print('{:3d}'.format(int(us[k].value())), end=" ")
print("")
# Tkinter による経路の描画
# draw(size, xs)
ハミルトン閉路のプログラムに制約条件を追加しただけなので、難しいところはないと思います。
●実行結果
実行結果は次のようになりました。
>>> import knight2 >>> import time >>> s = time.time(); knight2.solver(6); time.time() - s Status Optimal 0 7 34 15 24 9 33 16 1 8 31 14 6 35 32 25 10 23 17 26 19 2 13 30 20 5 28 11 22 3 27 18 21 4 29 12 0.23624825477600098 >>> s = time.time(); knight2.solver(8); time.time() - s Status Optimal 0 21 18 3 16 53 8 5 19 48 1 54 7 4 15 52 22 63 20 17 2 55 6 9 47 42 49 24 37 14 51 56 62 23 38 43 50 31 10 33 41 46 25 36 13 34 57 30 26 61 44 39 28 59 32 11 45 40 27 60 35 12 29 58 0.5663492679595947 >>> s = time.time(); knight2.solver(10); time.time() - s Status Optimal 0 81 92 89 2 83 30 41 4 43 93 96 1 82 91 88 3 44 29 40 80 99 90 95 68 31 84 87 42 5 97 94 79 50 85 48 45 6 39 28 78 51 98 67 32 69 86 17 46 7 55 66 77 60 49 36 47 8 27 38 52 59 54 35 70 33 18 37 16 9 65 56 73 76 61 22 15 12 19 26 74 53 58 63 34 71 24 21 10 13 57 64 75 72 23 62 11 14 25 20 66.20720839500427 >>> s = time.time(); knight2.solver(12); time.time() - s Status Optimal 0 121 76 139 2 141 74 59 70 57 84 63 77 138 1 142 75 60 71 56 73 62 69 86 122 143 120 137 140 3 8 61 58 85 64 83 135 78 123 116 9 80 55 72 65 82 87 68 124 119 136 79 4 7 10 81 20 67 16 13 111 134 117 6 115 54 21 66 11 14 19 88 118 125 112 53 92 5 94 25 22 17 12 15 133 110 127 114 95 24 91 102 97 26 89 18 126 113 132 109 52 93 96 23 90 29 98 27 131 128 49 46 43 40 103 106 101 36 33 30 48 45 130 51 108 105 38 41 34 31 28 99 129 50 47 44 39 42 107 104 37 100 35 32 16.84189200401306
- 実行環境 : Python 3.10.6, PuLP 2.7.0, Ubunts 22.04 LTS (WSL2, Windows 10), Intel Core i5-6200U 2.30GHz
MTZ 方式の制約を強化すると、実行速度はもう少し速くなるようです。制約の強化は「巡回セールスマン問題」で取り上げる予定です。
12 行 12 列盤の解の一例
●プログラムリスト
#
# knight2.py : 騎士の周遊
#
# Copyright (C) 2019-2023 Makoto Hiroi
#
import pulp
import tkinter as tk
# 描画
def draw(size, xs):
root = tk.Tk()
c0 = tk.Canvas(root, width = 40 * size + 20, height = 40 * size + 20)
c0.pack()
for y in range(size):
for x in range(size):
if (y + x) % 2 == 0:
color = 'white'
else:
color = 'black'
x1 = x * 40 + 10
y1 = y * 40 + 10
c0.create_rectangle(x1, y1, x1 + 40, y1 + 40, fill = color)
for i in range(size ** 2):
for j in range(size ** 2):
if xs[i][j].value() == 0: continue
y1 = (i // size) * 40 + 30
x1 = (i % size) * 40 + 30
y2 = (j // size) * 40 + 30
x2 = (j % size) * 40 + 30
c0.create_line(x1, y1, x2, y2, fill = 'cyan', width = 3.0)
root.mainloop()
# 隣接行列の生成
def make_matrix(size):
size2 = size * size
ws = [[0] * size2 for _ in range(size2)]
for y in range(size):
for x in range(size):
k = y * size + x
for dx, dy in [(-1, -2), (1, -2), (-2, -1), (2, -1),
(-2, 1), (2, 1), (-1, 2), (1, 2)]:
x1 = x + dx
y1 = y + dy
if 0 <= x1 < size and 0 <= y1 < size:
ws[k][y1 * size + x1] = 1
return ws
def solver(size):
size2 = size * size
ws = make_matrix(size) # 隣接行列
prob = pulp.LpProblem('knight-circle')
# 変数の生成
xs = []
for i in range(size2):
xs1 = [pulp.LpVariable('x{}_{}'.format(i, j), cat='Binary') for j in range(size2)]
xs.append(xs1)
us = [pulp.LpVariable('u{}'.format(i), cat='Integer', lowBound=0, upBound=size2-1) for i in range(size2)]
# 制約条件
for i in range(size2):
prob += pulp.lpSum([xs[j][i] for j in range(size2)]) == 1
prob += pulp.lpSum([xs[i][j] for j in range(size2)]) == 1
# 重み 0 の辺は選択しない
for i in range(size2):
for j in range(size2):
if ws[i][j] == 0: prob += xs[i][j] == 0
for i in range(size2):
for j in range(1, size2):
if i == j: continue
prob += us[i] + 1 - (size2 - 1) * (1 - xs[i][j]) <= us[j]
prob += us[0] == 0 # スタート
for i in range(1, size2): 1 <= us[i]
prob += xs[0][size + 2] == 1 # 最初に選ぶ辺を限定
prob += us[size + 2] == 1
prob += xs[size * 2 + 1][0] == 1 # 最後に選ぶ辺
prob += us[size * 2 + 1] == size2 - 1 # 最後の頂点
status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for i in range(size):
for j in range(size):
k = i * size + j
print('{:3d}'.format(int(us[k].value())), end=" ")
print("")
# Tkinter による経路の描画
# draw(size, xs)