Python3 Programming

●ビンパッキング問題

「ビンパッキング問題 (Bin Packing Problem)」は、大きさが異なる N 個の品物を大きさ B のビンに詰めるとき、最小のビンの本数と品物の詰め方を求める問題です。ビンは品物を入れる箱と考えてください。品物の大きさは B 以下で、ビンも品物を詰めるだけの本数は十分に用意されているものとします。

ビンパッキング問題は部分和問題と同じく NP 問題になるので、品物の数が多くなると最適解を求めるのが大変難しくなります。そこで今回は、ビンの本数をあらかじめ決めておいて、品物の詰め方を求めるプログラムを作ってみましょう。それでは問題です。

このような問題を解く場合、データ構造に「結合行列」を使うと簡単です。次の図を見てください。

   : a b c d e f g h
---+-----------------
 A : 1 0 1 0 0 0 0 0
 B : 0 0 0 1 0 0 1 0
 C : 0 1 0 0 1 0 0 0
 D : 0 0 0 0 0 1 0 1


    図 : 結合行列

a - h が品物、A - D が人とします。結合行列は人が運ぶ品物を 1 で、運ばない品物を 0 で表します。上図の場合、A が運ぶ品物は a, c で、D が運ぶ品物が f, h となります。制約の記述も簡単です。一人が運ぶ品物の重さの制限は lpDot() を使えば簡単ですね。もう一つ制約があって、品物を表す列には 1 が一つしかないことです。

この場合、lpSum() を使うと簡単に制約を表すことができます。たとえば、リスト xs の要素が 0 と 1 だけの場合、lpSum(xs) == 1 とすれば、xs には 1 が一つしかないという制約を表すことができます。

プログラムは次のようになります。

リスト : ビンパッキング問題

import pulp

# 人は４人
ps = ['A', 'B', 'C', 'D']

# 品物の重さ
ws = [3.3, 6.1, 5.8, 4.1, 5.0, 2.1, 6.0, 6.4]

prob = pulp.LpProblem('binpacking', sense = pulp.LpMaximize)

# 制約条件
vs = []
for p in ps:
    vs1 = [pulp.LpVariable('{}{}'.format(p, x), cat='Binary') for x in range(len(ws))]
    vs.append(vs1)
    prob += pulp.lpDot(ws, vs1) <= 11.0

for x in range(len(ws)):
    xs = [vs[i][x] for i in range(len(ps))]
    prob += pulp.lpSum(xs) == 1

print(prob)

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
print(ws)
for vs1 in vs:
    print([v.value() for v in vs1])

変数 vs に結合行列をセットします。一人につき変数 (Binary, 0 or 1) を 8 個生成してリストに格納し、それを vs に追加していきます。一人で 11.0 kg まで持つことができるので、制約条件は lpDot(ws, vs1) <= 11.0 となります。結合行列が完成したら列の制約条件を設定します。x 列目の変数をリストに格納して変数 xs にセットします。この合計値が 1 であればいいので、条件は lpSum(xs) == 1 となります。

それでは実行してみましょう。

binpacking:
MAXIMIZE
None
SUBJECT TO
_C1: 3.3 A0 + 6.1 A1 + 5.8 A2 + 4.1 A3 + 5 A4 + 2.1 A5 + 6 A6 + 6.4 A7 <= 11

_C2: 3.3 B0 + 6.1 B1 + 5.8 B2 + 4.1 B3 + 5 B4 + 2.1 B5 + 6 B6 + 6.4 B7 <= 11

_C3: 3.3 C0 + 6.1 C1 + 5.8 C2 + 4.1 C3 + 5 C4 + 2.1 C5 + 6 C6 + 6.4 C7 <= 11

_C4: 3.3 D0 + 6.1 D1 + 5.8 D2 + 4.1 D3 + 5 D4 + 2.1 D5 + 6 D6 + 6.4 D7 <= 11

_C5: A0 + B0 + C0 + D0 = 1

_C6: A1 + B1 + C1 + D1 = 1

_C7: A2 + B2 + C2 + D2 = 1

_C8: A3 + B3 + C3 + D3 = 1

_C9: A4 + B4 + C4 + D4 = 1

_C10: A5 + B5 + C5 + D5 = 1

_C11: A6 + B6 + C6 + D6 = 1

_C12: A7 + B7 + C7 + D7 = 1

VARIABLES
0 <= A0 <= 1 Integer
0 <= A1 <= 1 Integer
0 <= A2 <= 1 Integer

  ・・・省略・・・

0 <= D5 <= 1 Integer
0 <= D6 <= 1 Integer
0 <= D7 <= 1 Integer

Status Optimal
[3.3, 6.1, 5.8, 4.1, 5.0, 2.1, 6.0, 6.4]
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0]
[0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0]
[0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0]

A が 9.3 kg, B が 8.5 kg, C が 10.8 kg, D が 10.2 kg となりました。総重量 (38.8 kg) を 4 人で割ると 9.7 kg になりますが、これを上限値とすると、解はなくなります。4 人の場合、解があるのは 10.8 kg までです。興味のある方はいろいろ試してみてください。

●N Queens Problem

N Queens Problem は「8 クイーン」の拡張バージョンで、N 行 N 列の盤面に N 個のクイーンを互いの利き筋が重ならないように配置する問題です。クイーンは将棋の飛車と角をあわせた駒で、縦横斜めに任意に動くことができます。8 クイーンの解答例を示しましょう。

        図 : 8 クイーンの解答例

N Queens Problem の場合、行と列にはクイーンが必ず一つ存在し、斜め方向にはクイーンが一つ存在するラインと、クイーンが存在しないラインがあります。PuLP (CBC) で N Queens Problem を解くときは、これらを制約条件として使います。N 行 N 列の配列に変数 (Binary) をセットします。すると制約条件は、行方向の変数の和 == 1, 列方向の変数の和 == 1, 斜め方向の変数の和 <= 1 と表すことができます。

これをプログラムすると次のようになります。

リスト : N Queens Problem

import pulp

# 左上から右下へ
def make_line1(x, y, n):
    xs = []
    while x < n and y < n:
        xs.append((y, x))
        x += 1
        y += 1
    return xs

# 右上から左下へ
def make_line2(x, y, n):
    ys = []
    while x >= 0 and y < n:
        ys.append((y, x))
        x -= 1
        y += 1
    return ys

def solver(n):
    problem = pulp.LpProblem('sample')

    # 変数の生成
    vs = []
    for y in range(n):
        vs1 = [pulp.LpVariable('v_{}_{}'.format(y, x), cat='Binary') for x in range(n)]
        vs.append(vs1)

    # 行と列の制約条件
    for y in range(n):
        problem += pulp.lpSum([vs[y][x] for x in range(n)]) == 1
        problem += pulp.lpSum([vs[x][y] for x in range(n)]) == 1

    # 斜め方向の制約条件
    for z in range(n):
        problem += pulp.lpSum([vs[y][x] for y, x in make_line1(z, 0, n)]) <= 1
        problem += pulp.lpSum([vs[y][x] for y, x in make_line2(z, 0, n)]) <= 1
    for z in range(1, n):
        problem += pulp.lpSum([vs[y][x] for y, x in make_line1(0, z, n)]) <= 1
        problem += pulp.lpSum([vs[y][x] for y, x in make_line2(n - 1, z, n)]) <= 1
    # print(problem)

    status = problem.solve(pulp.PULP_CBC_CMD(msg=0))
    print("Status", pulp.LpStatus[status])
    for y in range(n):
        for x in range(n):
            if vs[y][x].value() == 1:
                print("Q", end=" ")
            else:
                print(".", end=" ")
        print("")

最初に変数 (Binary) を生成して配列 vs にセットします。行と列の制約条件は簡単です。斜め方向の場合、関数 make_line1 と make_line2 を使って座標を生成し、その位置にある変数を内包表記で取り出してリストに格納します。あとは、その和が 1 以下であることを設定すれば OK です。

まずは最初に、 4 * 4 盤で制約条件が正しく設定されているかチェックします。

>>> import nqueen
>>> nqueen.solver(4)
sample:
MINIMIZE
None
SUBJECT TO
_C1: v_0_0 + v_0_1 + v_0_2 + v_0_3 = 1

_C2: v_0_0 + v_1_0 + v_2_0 + v_3_0 = 1

_C3: v_1_0 + v_1_1 + v_1_2 + v_1_3 = 1

_C4: v_0_1 + v_1_1 + v_2_1 + v_3_1 = 1

_C5: v_2_0 + v_2_1 + v_2_2 + v_2_3 = 1

_C6: v_0_2 + v_1_2 + v_2_2 + v_3_2 = 1

_C7: v_3_0 + v_3_1 + v_3_2 + v_3_3 = 1

_C8: v_0_3 + v_1_3 + v_2_3 + v_3_3 = 1

_C9: v_0_0 + v_1_1 + v_2_2 + v_3_3 <= 1

_C10: v_0_0 <= 1

_C11: v_0_1 + v_1_2 + v_2_3 <= 1

_C12: v_0_1 + v_1_0 <= 1

_C13: v_0_2 + v_1_3 <= 1

_C14: v_0_2 + v_1_1 + v_2_0 <= 1

_C15: v_0_3 <= 1

_C16: v_0_3 + v_1_2 + v_2_1 + v_3_0 <= 1

_C17: v_1_0 + v_2_1 + v_3_2 <= 1

_C18: v_1_3 + v_2_2 + v_3_1 <= 1

_C19: v_2_0 + v_3_1 <= 1

_C20: v_2_3 + v_3_2 <= 1

_C21: v_3_0 <= 1

_C22: v_3_3 <= 1

VARIABLES
0 <= v_0_0 <= 1 Integer
0 <= v_0_1 <= 1 Integer
0 <= v_0_2 <= 1 Integer
0 <= v_0_3 <= 1 Integer
0 <= v_1_0 <= 1 Integer
0 <= v_1_1 <= 1 Integer
0 <= v_1_2 <= 1 Integer
0 <= v_1_3 <= 1 Integer
0 <= v_2_0 <= 1 Integer
0 <= v_2_1 <= 1 Integer
0 <= v_2_2 <= 1 Integer
0 <= v_2_3 <= 1 Integer
0 <= v_3_0 <= 1 Integer
0 <= v_3_1 <= 1 Integer
0 <= v_3_2 <= 1 Integer
0 <= v_3_3 <= 1 Integer

Status Optimal
. Q . .
. . . Q
Q . . .
. . Q .

行と列で制約条件が 8 本、斜め方向の制約条件が 7 * 2 = 14 本、制約条件は合計で 24 本になります。正常に動作していますね。

それでは盤面を大きくして試してみましょう。

>>> import time
>>> s = time.time(); nqueen.solver(8); time.time() - s
Status Optimal
. . . . . Q . .
. Q . . . . . .
. . . . . . Q .
Q . . . . . . .
. . Q . . . . .
. . . . Q . . .
. . . . . . . Q
. . . Q . . . .
0.018883705139160156
>>> s = time.time(); nqueen.solver(16); time.time() - s
Status Optimal
. . . . . Q . . . . . . . . . .
. . . . . . . . . . . Q . . . .
. . Q . . . . . . . . . . . . .
. . . . . . Q . . . . . . . . .
. . . . . . . . . . . . Q . . .
. Q . . . . . . . . . . . . . .
. . . . . . . . . . . . . Q . .
. . . . Q . . . . . . . . . . .
. . . . . . . . . . Q . . . . .
. . . . . . . . Q . . . . . . .
Q . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . Q
. . . Q . . . . . . . . . . . .
. . . . . . . . . . . . . . Q .
. . . . . . . . . Q . . . . . .
. . . . . . . Q . . . . . . . .
0.05153059959411621
>>> s = time.time(); nqueen.solver(32); time.time() - s
Status Optimal
. . . . . . . . . . . . . . . . . . . . . . Q . . . . . . . . .
. . . . . . . . . . . . . Q . . . . . . . . . . . . . . . . . .
. . . . . . . . Q . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . Q . . .
. . . . . . . . . . . . . . . Q . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q . .
. . . . . . . . . Q . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . Q . . . . . . . . . . .
. Q . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q .
. . . . . . . Q . . . . . . . . . . . . . . . . . . . . . . . .
. . Q . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . Q . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . Q . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . Q . . . . . . . . . .
. . . Q . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . Q . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . Q . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . Q . . . . .
. . . . . . . . . . . . . . Q . . . . . . . . . . . . . . . . .
. . . . . . Q . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q
. . . . . . . . . . . . . . . . . . . . . . . . . . . Q . . . .
Q . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . Q . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . Q . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . Q . . . . . .
. . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . Q . . . . . . . . . . . . .
. . . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . Q . . . . . . .
. . . . . . . . . . . . . . . . Q . . . . . . . . . . . . . . .
0.2751781940460205

• 実行環境 : Python 3.10.6, PuLP 2.7.0, Ubunts 22.04 LTS (WSL2, Windows 10), Intel Core i5-6200U 2.30GHz

32 * 32 盤でも 1 秒かからずに解くことができました。ちなみに、100 * 100 盤では 3.4 秒かかりました。CBC は非商用 (無料) ですが優秀なソルバーだと思いました。

●ナンバープレース

次は皆さんお馴染みのパズル「数独 (ナンバープレース)」の解法プログラムを作りましょう。ナンバープレースは 9×9 の盤を用いて、縦 9 列、横 9 列のそれぞれに 1 から 9 までの数字をひとつずつ入れます。また、太線で囲まれた 3×3 の枠内にも 1 から 9 までの数字をひとつずつ入れます。ただし、縦、横、枠の中で、同じ数字が重複して入ることはありません。

ナンバープレースの盤面 (9 行 9 列) を下図に示します。

  列 0 1 2   3 4 5   6 7 8
行 +-------+-------+-------+
 0 |       |       |       |
 1 | 枠 0  |   1   |   2   |
 2 |       |       |       |
   +-------+-------+-------+
 3 |       |       |       |
 4 |   3   |   4   |   5   |
 5 |       |       |       |
   +-------+-------+-------+
 6 |       |       |       |
 7 |   6   |   7   |   8   |
 8 |       |       |       |
   +-------+-------+-------+

   図 : 数独 (9 * 9) の盤面

この盤面を次のように二次元配列 (リストのリスト) で表すことにします。

リスト : 問題 (出典: 数独 - Wikipedia の問題例)

q00 = [
    [5, 3, 0,  0, 7, 0,  0, 0, 0],
    [6, 0, 0,  1, 9, 5,  0, 0, 0],
    [0, 9, 8,  0, 0, 0,  0, 6, 0],

    [8, 0, 0,  0, 6, 0,  0, 0, 3],
    [4, 0, 0,  8, 0, 3,  0, 0, 1],
    [7, 0, 0,  0, 2, 0,  0, 0, 6],

    [0, 6, 0,  0, 0, 0,  2, 8, 0],
    [0, 0, 0,  4, 1, 9,  0, 0, 5],
    [0, 0, 0,  0, 8, 0,  0, 7, 9]
]

各マスに一つの変数を対応させると、変数の値は 1 から 9 までになります。ところが、この方法では PuLP で解くことはできません。この場合、行と列と枠の 9 つの変数で異なる数字が一つずつ入ることになりますが、「変数の値がすべて異なる」という制約条件を表す方法が PuLP にはないからです。たとえば、SWI-Prolog の制約論理プログラミング clpfd には all_different という述語が用意されていて、「変数の値がすべて異なる」という制約条件を記述することが可能です。

そこで、各マスに 9 つの変数 (Binary) を用意して、変数と数字を対応させることにします。この場合、制約条件は次のように記述することができます。

1. 各マスに置くことができる数字は一つだけ
2. 1 から 9 までの各数字は行、列、枠で一つだけ

上記の条件に合う変数を取り出して、その合計値が 1 であれば制約条件を満たすことになります。これをプログラムすると次のようになります。

リスト : ナンバープレースの解法

import pulp

#
# 問題は省略
#

prob = pulp.LpProblem('numplace')

# 変数の生成
vs = []
for y in range(9):
    vs1 = []
    for x in range(9):
        vs2 = [pulp.LpVariable('v{}_{}_{}'.format(y, x, i), cat='Binary') for i in range(1, 10)]
        vs1.append(vs2)
    vs.append(vs1)

# 各マスの制約
for y in range(9):
    for x in range(9):
        prob += pulp.lpSum([vs[y][x][i] for i in range(9)]) == 1

# 行と列の制約
for j in range(9):
    for i in range(9):
        prob += pulp.lpSum([vs[j][n][i] for n in range(9)]) == 1
        prob += pulp.lpSum([vs[n][j][i] for n in range(9)]) == 1

# 枠の生成
bs = []
for y in [0, 3, 6]:
    for x in [0, 3, 6]:
        bs1 = [(y + i, x + j) for i in range(3) for j in range(3)]
        bs.append(bs1)

# 枠の制約
for b in bs:
    for i in range(9):
        prob += pulp.lpSum([vs[y][x][i] for y, x in b]) == 1

# 問題の設定
for y in range(9):
    for x in range(9):
        n = q00[y][x]
        if n == 0: continue
        prob += vs[y][x][n - 1] == 1

# print(prob)

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
print("Result")
for y in range(9):
    for x in range(9):
        print(sum([int(vs[y][x][n - 1].value() * n) for n in range(1, 10)]), end = " ")
    print("")

生成した変数 (Binary) は配列 vs に格納します。vs は三次元配列になることに注意してください。つまり、マス (y, x) の数字 n に対応する変数は vs[y][x][n - 1] になります。次に、マスと行列枠の制約条件を設定します。マスの制約条件は簡単ですね。行列枠の場合は同じ数字に対応する変数を取り出して、その合計値が 1 であることを設定します。最後に問題を読み込み、該当するマスの数字を 1 に設定するだけです。

それでは実行してみましょう。

Status Optimal
Result
5 3 4 6 7 8 9 1 2
6 7 2 1 9 5 3 4 8
1 9 8 3 4 2 5 6 7
8 5 9 7 6 1 4 2 3
4 2 6 8 5 3 7 9 1
7 1 3 9 2 4 8 5 6
9 6 1 5 3 7 2 8 4
2 8 7 4 1 9 6 3 5
3 4 5 2 8 6 1 7 9

一瞬で答えを求めることができました。

●覆面算

計算式の数字を文字や記号に置き換えて、それを元の数字に戻すパズルを「覆面算」といいます。異なる文字は異なる数字を表し、同じ文字は同じ数字を表します。使用する数字は 0 から 9 までで、最上位の桁に 0 を入れることはできません。

     S E N D
 +   M O R E
-------------
   M O N E Y

  図 : 覆面算

問題はデュードニーが 1924 年に発表したもので、覆面算の古典といわれる有名なパズルです。

それではプログラムを作りましょう。文字は s, e, n, d, m, o, r, y の 8 つあり、そこに異なる数字 (0 - 9) を割り当てます。PuLP には変数の値がすべて異なる (たとえば clpfd の all_different) という制約条件は記述できないので、各文字に 10 個の変数 (Binary) を用意して、それと数字を対応させることにします。すると、制約条件は次のように記述することができます。

1. 記号に対応する変数の総和は 1 になる
2. 数字 0 - 9 に対応する変数の総和は 1 以下になる
3. 式 send + more == money が成り立つ
4. send + more = money は足し算なので、m は 1 になる

今回の問題は 10 個数字を 8 個の文字に割り当てるので、使用しない数字もあります。このため、2 の条件が <= 1 になります。== 1 とすると解くことができません。ご注意くださいませ。

プログラムは次のようになります。

リスト : 覆面算

import pulp

prob = pulp.LpProblem('hukumenn')

# 数字
nums = list(range(10))

# 変数の生成
s = [pulp.LpVariable('s{}'.format(i), cat='Binary') for i in range(10)]
e = [pulp.LpVariable('e{}'.format(i), cat='Binary') for i in range(10)]
n = [pulp.LpVariable('n{}'.format(i), cat='Binary') for i in range(10)]
d = [pulp.LpVariable('d{}'.format(i), cat='Binary') for i in range(10)]
m = [pulp.LpVariable('m{}'.format(i), cat='Binary') for i in range(10)]
o = [pulp.LpVariable('o{}'.format(i), cat='Binary') for i in range(10)]
r = [pulp.LpVariable('r{}'.format(i), cat='Binary') for i in range(10)]
y = [pulp.LpVariable('y{}'.format(i), cat='Binary') for i in range(10)]

# 制約条件
for x in [s, e, n, d, m, o, r, y]:
    prob += pulp.lpSum(x) == 1
for x in range(10):
    prob += pulp.lpSum([v[x] for v in [s, e, n, d, m, o, r, y]]) <= 1

send = 1000 * pulp.lpDot(s, nums) + 100 * pulp.lpDot(e, nums) \
     + 10 * pulp.lpDot(n, nums) + pulp.lpDot(d, nums)
more = 1000 * pulp.lpDot(m, nums) + 100 * pulp.lpDot(o, nums) \
     + 10 * pulp.lpDot(r, nums) + pulp.lpDot(e, nums)
money = 10000 * pulp.lpDot(m, nums) + 1000 * pulp.lpDot(o, nums) + 100 * pulp.lpDot(n, nums) \
      + 10 * pulp.lpDot(e, nums) + pulp.lpDot(y, nums)

prob += send + more == money
prob += m[1] == 1

status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
print("Status", pulp.LpStatus[status])
for x in [s, e, n, d, m, o, r, y]:
    print(sum([i.value() * j for i, j in zip(x, nums)]), end=" ")
print("")

制約条件をそのままプログラムしただけなので、とくに難しいところはないと思います。それでは実行してみましょう。

Status Optimal
9.0 5.0 6.0 7.0 1.0 0.0 8.0 2.0

答えは 9567 + 1085 = 10652 となりました。実行時間は M.Hiroi の環境で 0.37 秒でした。ところで、制約条件 prob += m[1] == 1 のかわりに、prob += s[0] == 0, prob += m[0] == 0 としても解くことができます。興味のある方は試してみてください。

●魔方陣

次は、魔方陣を解くプログラムを作ってみましょう。

          図 : 魔方陣

上図は 3 行 3 列の魔方陣ですが、N 行 N 列の魔方陣を解くプログラムを作ります。ただし、N が大きくなると現実的な時間では解けないと思います。ご注意くださいませ。

魔方陣 - Wikipedia によると、n * n の正方形の方陣に 1 から n² の数字を配置するとき、一列の和 w は次式で計算できるそうです。

w = (1 + 2 + 3 + ... + n2) / n = n * (n2 + 1) / 2

今回も各マスに n² 個の変数 (Binary) を用意して、制約条件を記述することにします。

1. 各マスに対応する変数の総和は 1 になる
2. 数字 1 - n² に対応する変数の総和は 1 になる
3. 行、列、斜めの数字の和が w になる

プログラムは次のようになります。

リスト : 魔方陣 (magic.py)

import pulp
from itertools import product

def solver(size):
    # 数字
    nums = list(range(1, size * size + 1))
    w = size * (size * size + 1) // 2

    prob = pulp.LpProblem('magic-square')
    # 変数
    board = []
    for y in range(size):
        bs1 = []
        for x in range(size):
            bs = [pulp.LpVariable('b{}_{}_{}'.format(y, x, n), cat='Binary') for n in nums]
            bs1.append(bs)
        board.append(bs1)

    # 制約条件
    # 各マスの数字は一つ
    for x, y in product(range(size), range(size)):
        prob += pulp.lpSum([board[y][x][n] for n in range(len(nums))]) == 1
    # 数字は盤面全体で一つだけ
    for n in range(len(nums)):
        prob += pulp.lpSum([board[y][x][n] for x, y in product(range(size), range(size))]) == 1
    # 縦
    for y in range(size):
        bs = [board[y][x] for x in range(size)]
        prob += pulp.lpSum([pulp.lpDot(b, nums) for b in bs]) == w
    # 横
    for x in range(size):
        bs = [board[y][x] for y in range(size)]
        prob += pulp.lpSum([pulp.lpDot(b, nums) for b in bs]) == w
    # 斜め
    bs = [board[y][x] for x, y in zip(range(size), range(size))]
    prob += pulp.lpSum([pulp.lpDot(b, nums) for b in bs]) == w
    bs = [board[y][x] for x, y in zip(range(size), range(size - 1, -1, -1))]
    prob += pulp.lpSum([pulp.lpDot(b, nums) for b in bs]) == w

    #print(prob)
    status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
    print("Status", pulp.LpStatus[status])
    for y in range(size):
        for x in range(size):
            print(int(sum([a.value() * n for a, n in zip(board[y][x], nums)])), end = " ")
        print("")

プログラムは簡単なので説明は割愛いたします。それでは実際に試してみましょう。

>>> import magic
>>> import time
>>> s = time.time(); magic.solver(3); time.time() - s
Status Optimal
4 3 8
9 5 1
2 7 6
0.037772178649902344
>>> s = time.time(); magic.solver(4); time.time() - s
Status Optimal
8 10 15 1
4 11 14 5
9 6 3 16
13 7 2 12
4.930948257446289

• 実行環境 : Python 3.10.6, PuLP 2.7.0, Ubunts 22.04 LTS (WSL2, Windows 10), Intel Core i5-6200U 2.30GHz

5 行 5 列の魔方陣は時間がかかるようで、途中であきらめました。そこで、対称解のチェックを入れて解を制限してみることにします。対称解のチェックは面倒だと思われる方もいるでしょう。ところが、下図のように四隅の大小関係を利用すると簡単です。

      図 : 対称解のチェック

魔方陣の場合、回転解が 4 種類あって、鏡像解が 2 種類あります。四隅の大小関係をチェックすることで、これらの対称解を排除することができます。プログラムは制約条件を追加するだけです。

リスト :  対称解のチェック

    # A B C   A < C < G
    # D E F   A < I
    # G H I
    prob += pulp.lpSum(pulp.lpDot(board[0][0], nums)) <= pulp.lpSum(pulp.lpDot(board[0][size - 1], nums))
    prob += pulp.lpSum(pulp.lpDot(board[0][size - 1], nums)) <= pulp.lpSum(pulp.lpDot(board[size - 1][0], nums))
    prob += pulp.lpSum(pulp.lpDot(board[0][0], nums)) <= pulp.lpSum(pulp.lpDot(board[size - 1][size - 1], nums))

実行結果は次のようになりました。

>>> import magic1
>>> import time
>>> s = time.time(); magic1.solver(3); time.time() - s
Status Optimal
2 9 4
7 5 3
6 1 8
0.09526181221008301
>>> s = time.time(); magic1.solver(4); time.time() - s
Status Optimal
7 4 14 9
11 16 2 5
6 13 3 12
10 1 15 8
2.7991225719451904

制約条件を追加しましたが、5 行 5 列の魔方陣は解くことができませんでした。参考文献 には、n 行 n 列の魔方陣で n が奇数の場合の簡単な作成プログラムが紹介されています。興味のある方はいろいろ試してみてください。

●最短経路問題

次は下図に示す経路で、A から G までの最短経路 (最小コストの経路) を求めてみましょう。

      図 : 経路図

 表 : 辺のコスト

   辺  : コスト
 ------+--------
 A - B :    1
 A - C :    6
 B - C :    4
 B - D :    2
 C - E :    7
 D - E :    5
 D - F :    3
 E - G :    8

なお、最短経路問題は「ダイクストラのアルゴリズム」を使って高速に解くことができます。興味のある方は拙作のページ Algorithms with Python 欲張り法 をお読みください。

最短経路問題を最適化ソルバーで解く場合、辺を変数 (Binary) に対応させると簡単です。辺を選ぶときは変数の値を 1 とし、選ばない場合は 0 とします。頂点 i から頂点 j に行く辺を X_ij, その辺のコストを W_ij とすると、目的関数と制約条件は次の式で表すことができます。

目的関数:

Σ Σ Wij * Xij
i j

制約条件:

1. Σ Xai == 1, 頂点 a から出る辺は一つだけ
   i

2. Σ Xig == 1, 頂点 g に入る辺は一つだけ
   i

3. Σ Xij == Σ Xjk, 頂点 j に入る辺の数と j から出る辺の数は等しい (j != a, g)
   i        k

スタートとゴール以外の頂点の場合、辺がつながっていないときは入力が 0 で出力が 0 になります。辺がつながっているときは入力が 1 で出力が 1 になります。これを制約条件 3 で表しています。

プログラムは隣接行列を使うと簡単です。次のリストを見てください。

リスト : 最短経路問題

import pulp

#
# 辺の重み
#
#      A  B  C  D  E  F  G
ws = [[0, 1, 6, 0, 0, 0, 0],  # A
      [1, 0, 4, 2, 0, 0, 0],  # B
      [6, 4, 0, 0, 7, 0, 0],  # C
      [0, 2, 0, 0, 5, 3, 0],  # D
      [0, 0, 7, 5, 0, 0, 8],  # E
      [0, 0, 0, 3, 0, 0, 0],  # F
      [0, 0, 0, 0, 8, 0, 0],  # G
]

node = ['A', 'B', 'C', 'D', 'E', 'F', 'G']
size  = len(node)   # 頂点の数
start = 0           # スタート (A)
goal  = 5           # ゴール (G)

prob = pulp.LpProblem('keiro')

# 変数の生成
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([pulp.lpDot(w, x) for w, x in zip(ws, xs)])

# 制約条件
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)]) == pulp.lpSum([xs[i][j] for j in range(size)])

# 重み 0 の辺は選択しない
for i in range(size):
    for j in range(size):
        if ws[i][j] == 0: prob += xs[i][j] == 0

# print(prob)

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("z", prob.objective.value())

辺の重み (コスト) は配列 ws に格納します。ws[i][j] が W_ij を表します。0 は辺がないことを表します。変数 (Binary) は配列 xs にセットします。xs[i][j] が変数 X_ij を表します。目的関数と制約条件は式をそのままプログラムするだけです。あとは ws[i][j] が 0 であれば、制約条件に xs[i][j] == 0 を追加します。

それでは実行してみましょう。

Status Optimal
A -> B
B -> D
D -> E
E -> G
z 16.0

最短経路は A -> B -> D -> E -> G で、コストは 16 になりました。頂点の数が一番少なくなる経路は A - C - E - G ですが、これだとコストが 21 になってしまいます。ちなみに、ゴールを F にすると結果は次のようになります。

Status Optimal
A -> B
B -> D
D -> F
z 6.0

A -> B -> D -> F が最短経路でコストは 6 になります。これは頂点数が一番少ない経路と一致します。