「ビンパッキング問題 (Bin Packing Problem)」は、大きさが異なる N 個の品物を大きさ B のビンに詰めるとき、最小のビンの本数と品物の詰め方を求める問題です。ビンは品物を入れる箱と考えてください。品物の大きさは B 以下で、ビンも品物を詰めるだけの本数は十分に用意されているものとします。
ビンパッキング問題は部分和問題と同じく NP 問題になるので、品物の数が多くなると最適解を求めるのが大変難しくなります。そこで今回は、ビンの本数をあらかじめ決めておいて、品物の詰め方を求めるプログラムを作ってみましょう。それでは問題です。
4 人で 8 個の荷物を運びます。各荷物の重さは 3.3 kg, 6.1 kg, 5.8 kg, 4.1 kg, 5.0 kg, 2.1 kg, 6.0 kg, 6.4 kg です。各自の運ぶ荷物の重さの合計が 11 kg 以下になるように荷物を割り当てることはできるでしょうか。できる場合はその割り当て方を求めてください。
出典 : Coprisによる制約プログラミング入門, (田村直之さん)
このような問題を解く場合、データ構造に「結合行列」を使うと簡単です。次の図を見てください。
: 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 となります。制約の記述も簡単です。一人が運ぶ品物の重さの制限は sum() を使えば簡単ですね。もう一つ制約があって、品物を表す列には 1 が一つしかないことです。
この場合も sum() を使うと簡単に制約を表すことができます。たとえば、二次元配列 vs の要素が 0 と 1 だけの場合、sum(vs[j, i] for j = 1:4) == 1 とすれば、vs の i 列目には 1 が一つしかないという制約を表すことができます。
プログラムは次のようになります。
リスト : ビンパッキング問題
using JuMP, HiGHS
# 人は４人
ps = ['A', 'B', 'C', 'D']
# 品物の重さ
ws = [3.3, 6.1, 5.8, 4.1, 5.0, 2.1, 6.0, 6.4]
model = Model(HiGHS.Optimizer)
# 変数
@variable(model, vs[1:4, 1:8], Bin)
# 制約条件
@constraint(model, [i = 1:4], sum(vs[i, j] * ws[j] for j = 1:8) <= 11)
@constraint(model, [i = 1:8], sum(vs[j, i] for j = 1:4) == 1)
# 実行
println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果: ", termination_status(model))
println(value.(vs))
println(value.(vs) * ws)
二次元配列 vs に結合行列をセットします。変数の型は Bin (0 or 1) を指定します。一人で 11.0 kg まで持つことができるので、制約条件は sum(vs[i, j] * ws[j] for j = 1:8) <= 11 となります。列の制約条件も簡単で、i 番目の列は sum(vs[j, i] for j = 1:4) == 1 となります。
それでは実行してみましょう。
$ julia binpack.jl Feasibility Subject to vs[1,1] + vs[2,1] + vs[3,1] + vs[4,1] = 1 vs[1,2] + vs[2,2] + vs[3,2] + vs[4,2] = 1 vs[1,3] + vs[2,3] + vs[3,3] + vs[4,3] = 1 vs[1,4] + vs[2,4] + vs[3,4] + vs[4,4] = 1 vs[1,5] + vs[2,5] + vs[3,5] + vs[4,5] = 1 vs[1,6] + vs[2,6] + vs[3,6] + vs[4,6] = 1 vs[1,7] + vs[2,7] + vs[3,7] + vs[4,7] = 1 vs[1,8] + vs[2,8] + vs[3,8] + vs[4,8] = 1 3.3 vs[1,1] + 6.1 vs[1,2] + 5.8 vs[1,3] + 4.1 vs[1,4] + 5 vs[1,5] + 2.1 vs[1,6] + 6 vs[1,7] + 6.4 vs[1,8] ≤ 11 3.3 vs[2,1] + 6.1 vs[2,2] + 5.8 vs[2,3] + 4.1 vs[2,4] + 5 vs[2,5] + 2.1 vs[2,6] + 6 vs[2,7] + 6.4 vs[2,8] ≤ 11 3.3 vs[3,1] + 6.1 vs[3,2] + 5.8 vs[3,3] + 4.1 vs[3,4] + 5 vs[3,5] + 2.1 vs[3,6] + 6 vs[3,7] + 6.4 vs[3,8] ≤ 11 3.3 vs[4,1] + 6.1 vs[4,2] + 5.8 vs[4,3] + 4.1 vs[4,4] + 5 vs[4,5] + 2.1 vs[4,6] + 6 vs[4,7] + 6.4 vs[4,8] ≤ 11 vs[1,1] binary ... 略 ... vs[4,8] binary 結果: OPTIMAL [-0.0 1.0 -0.0 -0.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 -0.0 -0.0 1.0 -0.0 -0.0 -0.0 1.0; 1.0 -0.0 -0.0 -0.0 -0.0 -0.0 1.0 -0.0] [8.2, 10.8, 10.5, 9.3]
A が 8.2 kg, B が 10.8 kg, C が 10.5 kg, D が 9.3 kg となりました。総重量 (38.8 kg) を 4 人で割ると 9.7 kg になりますが、これを上限値とすると、解はなくなります。4 人の場合、解があるのは 10.8 kg までです。興味のある方はいろいろ試してみてください。
N Queens Problem は「8 クイーン」の拡張バージョンで、N 行 N 列の盤面に N 個のクイーンを互いの利き筋が重ならないように配置する問題です。クイーンは将棋の飛車と角をあわせた駒で、縦横斜めに任意に動くことができます。8 クイーンの解答例を示しましょう。
列
1 2 3 4 5 6 7 8
*-----------------*
1 | Q . . . . . . . |
2 | . . . . Q . . . |
3 | . . . . . . . Q |
行 4 | . . . . . Q . . |
5 | . . Q . . . . . |
6 | . . . . . . Q . |
7 | . Q . . . . . . |
8 | . . . Q . . . . |
*-----------------*
図 : 8 クイーンの解答例
N Queens Problem の場合、行と列にはクイーンが必ず一つ存在し、斜め方向にはクイーンが一つ存在するラインと、クイーンが存在しないラインがあります。JuMP で N Queens Problem を解くときは、これらを制約条件として使います。N 行 N 列の配列に変数 (Bin) をセットします。すると制約条件は、行方向の変数の和 == 1, 列方向の変数の和 == 1, 斜め方向の変数の和 <= 1 と表すことができます。
これをプログラムすると次のようになります。
リスト : N Queens Problem
using JuMP, HiGHS
# 盤面の表示
function print_board(xs, n)
for x in 1 : n
for y in 1 : n
if xs[x, y] == 1
print("Q ")
else
print(". ")
end
end
println("")
end
end
function writeline!(d, n, p)
if haskey(d, n)
push!(d[n], p)
else
d[n] = [p]
end
end
# 斜めの利き筋を生成する
function diagline(n)
dr = Dict()
dl = Dict()
for (x, y) = Iterators.product(1:n, 1:n)
writeline!(dr, x + y, (x, y))
writeline!(dl, x - y, (x, y))
end
(dr, dl)
end
function solver(n)
model = Model(HiGHS.Optimizer)
# 変数
@variable(model, vs[1:n, 1:n], Bin)
# 縦横の制約
for i = 1 : n
@constraint(model, sum(vs[i, j] for j = 1 : n) == 1)
@constraint(model, sum(vs[j, i] for j = 1 : n) == 1)
end
# 斜めの制約
(dr, dl) = diagline(n)
for xs = values(dr)
@constraint(model, sum(vs[x, y] for (x, y) in xs) <= 1)
end
for xs = values(dl)
@constraint(model, sum(vs[x, y] for (x, y) in xs) <= 1)
end
# 実行
println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果: ", termination_status(model))
print_board(value.(vs), n)
end
最初に変数 (Bin) を格納した二次元配列 vs を作ります。行と列の制約条件は簡単です。斜め方向の場合、関数 diagline() を使って座標を生成します。座標を (x, y) とすると、x - y が同じ値であれば、左上から右下への方向で同じ斜めラインに並んでいます。x + y が同じ値であれば、右上から左下への方向で同じ斜めラインに並んでいます。この関係を使って同じ斜めラインの座標を格納したベクタを生成します。diagline() は斜めラインを格納した辞書を返すので、辞書からラインを一つずつ取り出し、@constraint で制約条件 (和が 1 以下であること) を設定するだけです。
まずは最初に、 4 * 4 盤で制約条件が正しく設定されているかチェックします。
julia> include("nqueen.jl")
solver (generic function with 1 method)
julia> solver(4)
Feasibility
Subject to
vs[1,1] + vs[1,2] + vs[1,3] + vs[1,4] = 1
vs[1,1] + vs[2,1] + vs[3,1] + vs[4,1] = 1
vs[2,1] + vs[2,2] + vs[2,3] + vs[2,4] = 1
vs[1,2] + vs[2,2] + vs[3,2] + vs[4,2] = 1
vs[3,1] + vs[3,2] + vs[3,3] + vs[3,4] = 1
vs[1,3] + vs[2,3] + vs[3,3] + vs[4,3] = 1
vs[4,1] + vs[4,2] + vs[4,3] + vs[4,4] = 1
vs[1,4] + vs[2,4] + vs[3,4] + vs[4,4] = 1
vs[4,1] + vs[3,2] + vs[2,3] + vs[1,4] ≤ 1
vs[3,1] + vs[2,2] + vs[1,3] ≤ 1
vs[4,2] + vs[3,3] + vs[2,4] ≤ 1
vs[4,3] + vs[3,4] ≤ 1
vs[1,1] ≤ 1
vs[4,4] ≤ 1
vs[2,1] + vs[1,2] ≤ 1
vs[1,1] + vs[2,2] + vs[3,3] + vs[4,4] ≤ 1
vs[3,1] + vs[4,2] ≤ 1
vs[1,2] + vs[2,3] + vs[3,4] ≤ 1
vs[1,4] ≤ 1
vs[1,3] + vs[2,4] ≤ 1
vs[4,1] ≤ 1
vs[2,1] + vs[3,2] + vs[4,3] ≤ 1
vs[1,1] binary
vs[2,1] binary
vs[3,1] binary
vs[4,1] binary
vs[1,2] binary
vs[2,2] binary
vs[3,2] binary
vs[4,2] binary
vs[1,3] binary
vs[2,3] binary
vs[3,3] binary
vs[4,3] binary
vs[1,4] binary
vs[2,4] binary
vs[3,4] binary
vs[4,4] binary
結果: OPTIMAL
. Q . .
. . . Q
Q . . .
. . Q .
行と列で制約条件が 8 本、斜め方向の制約条件が 7 * 2 = 14 本、制約条件は合計で 22 本になります。正常に動作していますね。
それでは盤面を大きくして試してみましょう。
julia> @time solver(8) 結果: OPTIMAL . . Q . . . . . . . . . . . Q . . Q . . . . . . . . . . . . . Q . . . . Q . . . Q . . . . . . . . . . Q . . . . . . . . . Q . . 0.017528 seconds (7.20 k allocations: 302.750 KiB, 24.22% compilation time) julia> @time solver(16) 結果: OPTIMAL . Q . . . . . . . . . . . . . . . . . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . Q . . . . . . . Q . . . . . . . . . . . . . . . . . Q . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . . Q . . . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . . . . . . . Q . Q . . . . . . . . . . . . . . . . . . . . . . . . . . . . Q . . . . . . . Q . . . . . . . . . . . . . . . . . . . . . . Q . . . . . Q . . . . . . . . . . . . . . . . . . . Q . . . . . . . . . . . . . . . . . . . Q . . . . . 0.018234 seconds (16.63 k allocations: 886.008 KiB) julia> @time solver(32) 結果: 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.055360 seconds (51.53 k allocations: 2.894 MiB)
32 * 32 盤でも 0.1 秒かからずに解くことができました。ちなみに、100 * 100 盤でも 1 秒かからずに解くことができます。HiGHS は非商用 (無料) ですが優秀なソルバーだと思いました。
次は皆さんお馴染みのパズル「数独 (ナンバープレース)」の解法プログラムを作りましょう。ナンバープレースは 9×9 の盤を用いて、縦 9 列、横 9 列のそれぞれに 1 から 9 までの数字をひとつずつ入れます。また、太線で囲まれた 3×3 の枠内にも 1 から 9 までの数字をひとつずつ入れます。ただし、縦、横、枠の中で、同じ数字が重複して入ることはありません。
ナンバープレースの盤面 (9 行 9 列) を下図に示します。
列 1 2 3 4 5 6 7 8 9 行 +-------+-------+-------+ 1 | | | | 2 | 枠 1 | 2 | 3 | 3 | | | | +-------+-------+-------+ 4 | | | | 5 | 4 | 5 | 6 | 6 | | | | +-------+-------+-------+ 7 | | | | 8 | 7 | 8 | 9 | 9 | | | | +-------+-------+-------+ 図 : 数独 (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 までになります。ところが JuMP の場合、この方法では簡単に解くことができません。この場合、行と列と枠の 9 つの変数で異なる数字が一つずつ入ることになりますが、「変数の値がすべて異なる」という制約条件を表す簡単な方法が JuMP には用意されていないからです。たとえば、SWI-Prolog の制約論理プログラミング clpfd には all_different という述語が用意されていて、「変数の値がすべて異なる」という制約条件を簡単に記述することが可能です。
そこで、各マスに 9 つの変数 (Bin) を用意して、変数と数字を対応させることにします。この場合、制約条件は次のように記述することができます。
これをプログラムすると次のようになります。
リスト : ナンバープレースの解法
using JuMP, HiGHS
#
# 問題は省略
#
model = Model(HiGHS.Optimizer)
# 変数の定義
@variable(model, vs[1:9, 1:9, 1:9], Bin)
# マスの制約 (各マスに数字は一つ)
for x = 1 : 9
for y = 1 : 9
@constraint(model, sum(vs[x, y, n] for n = 1 : 9) == 1)
end
end
# 行と列の制約 (行に同じ数字は一つ, 列に同じ数字は一つ)
for n = 1 : 9
for x = 1 : 9
@constraint(model, sum(vs[x, y, n] for y = 1 : 9) == 1)
@constraint(model, sum(vs[y, x, n] for y = 1 : 9) == 1)
end
end
# 枠の制約 (枠に同じ数字は一つ)
for x = [1, 4, 7]
for y = [1, 4, 7]
zs = [(x1, y1) for x1 = x:x+2 for y1 = y:y+2]
for n = 1 : 9
@constraint(model, sum(vs[x1, y1, n] for (x1, y1) = zs) == 1)
end
end
end
# 問題の読み込み
for x = 1 : 9
for y = 1 : 9
n = q00[x, y]
if n != 0
@constraint(model, vs[x, y, n] == 1)
end
end
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果: ", termination_status(model))
for x = 1 : 9
for y = 1 : 9
print(floor(Int64, sum(value(vs[x, y, n]) * n for n = 1 : 9)), " ")
end
println("")
end
生成した変数 (Bin) は配列 vs に格納します。vs は三次元配列になることに注意してください。つまり、マス (x, y) の数字 n に対応する変数は vs[x, y, n] になります。次に、マスと行列枠の制約条件を設定します。マスの制約条件は簡単ですね。行列枠の場合は同じ数字に対応する変数を取り出して、その合計値が 1 であることを設定します。最後に問題を読み込み、該当するマスの数字を 1 に設定するだけです。
それでは実行してみましょう。
$ julia numplace.jl 結果: OPTIMAL 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) を割り当てます。ナンバープレースのように、各文字に 10 個の変数 (Bin) を用意して、それと数字を対応させる方法でも解くことができますが、文字に対応する変数は Int 型にしたほうが、わかりやすいプログラムになります。そこで、all different の制約は変数 (Bin) を格納した二次元配列を使ってを記述することにします。
最初に文字に対応する変数を定義しましょう。次のリストを見てください。
リスト : 変数の定義 # 変数 @variable(model, 1 <= s <= 9, Int) @variable(model, 0 <= e <= 9, Int) @variable(model, 0 <= n <= 9, Int) @variable(model, 0 <= d <= 9, Int) @variable(model, 1 <= m <= 9, Int) @variable(model, 0 <= o <= 9, Int) @variable(model, 0 <= r <= 9, Int) @variable(model, 0 <= y <= 9, Int) vs = [s, e, n, d, m, o, r, y] size = length(vs) # 式の定義 @expression(model, send, 1000s + 100e + 10n + d) @expression(model, more, 1000m + 100o + 10r + e) @expression(model, money, 10000m + 1000o + 100n + 10e + y) # 式の制約 @constraint(model, send + more == money)
変数 s, e, n, d, m, o, r, y を @variable で定義します。変数の範囲は 0 以上 9 以下ですが、s と m は 0 にならないので、下限値を 1 に設定します。これらの変数はベクタに格納して vs にセットします。次に、@expression で send, more, money の数式を定義します。すると式の制約は send + more == money と定義することができます。
次に、all different を表す制約を定義しましょう。次のリストを見てください。
リスト : all different を表す制約条件 @variable(model, ns[1:size, 1:10], Bin) @constraint(model, [v = 1:size], sum(ns[v, n] for n = 1:10) == 1) # 変数には必ず数字が一つ割り当てられる @constraint(model, [n = 1:10], sum(ns[v, n] for v = 1:size) <= 1) # 異なる数字を割り当てる # 変数と数字を結びつける @constraint(model, [v = 1:size], sum(ns[v, n+1] * n for n = 0:9) == vs[v])
二次元配列 ns に Bin 型変数を定義します。行 (v = 1:size) が変数 vs[v] に対応し、列 n が数字 0 - 9 に対応します。すると、all different は次の条件で表すことができます。
今回の問題は 10 個数字を 8 個の文字に割り当てるので、使用しない数字もあります。このため、2 の条件が <= 1 になります。== 1 とすると解くことができません。ご注意くださいませ。最後に、変数 vs と ns を関連付けます。ns[v, :] が表す数字は sum(ns[v, n+1] * n for n = 0:9) で求めることができるので、この値が vs[v] と等しくなればいいわけです。
あとはプログラムを実行するだけです。
$ julia hukumen.jl 結果: OPTIMAL 9.0 5.0 6.0 7.0 1.0 0.0 8.0 2.0
答えは 9567 + 1085 = 10652 となりました。
リスト : 覆面算
using JuMP, HiGHS
model = Model(HiGHS.Optimizer)
# 変数
@variable(model, 1 <= s <= 9, Int)
@variable(model, 0 <= e <= 9, Int)
@variable(model, 0 <= n <= 9, Int)
@variable(model, 0 <= d <= 9, Int)
@variable(model, 1 <= m <= 9, Int)
@variable(model, 0 <= o <= 9, Int)
@variable(model, 0 <= r <= 9, Int)
@variable(model, 0 <= y <= 9, Int)
vs = [s, e, n, d, m, o, r, y]
size = length(vs)
# 式の定義
@expression(model, send, 1000s + 100e + 10n + d)
@expression(model, more, 1000m + 100o + 10r + e)
@expression(model, money, 10000m + 1000o + 100n + 10e + y)
# 式の制約
@constraint(model, send + more == money)
# all different を表す制約
@variable(model, ns[1:size, 1:10], Bin)
@constraint(model, [v = 1:size], sum(ns[v, n] for n = 1:10) == 1) # 変数には必ず数字が一つ割り当てられる
@constraint(model, [n = 1:10], sum(ns[v, n] for v = 1:size) <= 1) # 異なる数字を割り当てる
# 変数と数字を結びつける
@constraint(model, [v = 1:size], sum(ns[v, n+1] * n for n = 0:9) == vs[v])
# 実行
#println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
for x in vs
print(sum(value(x)), " ")
end
次は、魔方陣を解くプログラムを作ってみましょう。
┌─┬─┬─┐ 式
│Ａ│Ｂ│Ｃ│ A + B + C = N, A + E + I = N
├─┼─┼─┤ D + E + F = N, C + E + G = N
│Ｄ│Ｅ│Ｆ│ G + H + I = N
├─┼─┼─┤ A + D + G = N
│Ｇ│Ｈ│Ｉ│ B + E + H = N
└─┴─┴─┘ C + F + I = N
図 : 魔方陣
上図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。縦横斜めの合計が等しくなるように数字を配置してください。
上図は 3 行 3 列の魔方陣ですが、N 行 N 列の魔方陣を解くプログラムを作ります。ただし、N が大きくなると現実的な時間では解けないと思います。ご注意くださいませ。
魔方陣 - Wikipedia によると、n * n の正方形の方陣に 1 から \(n^2\) の数字を配置するとき、一列の和 w は次式で計算できるそうです。
今回も覆面算と同様に、盤面の数字は Int 型変数で表し、二次元配列を使って all different の制約を表すことにします。プログラムは次のようになります。
リスト : 魔方陣 (magic.jl)
using JuMP, HiGHS
function solver(size)
model = Model(HiGHS.Optimizer)
# 数字
size2 = size * size
nums = 1 : size2
w = div(size * (size2 + 1), 2)
# 変数
@variable(model, 1 <= vs[1:size, 1:size] <= size2, Int)
# 魔法陣の制約
@constraint(model, [x = 1:size], sum(vs[x, y] for y = 1:size) == w) # 横
@constraint(model, [y = 1:size], sum(vs[x, y] for x = 1:size) == w) # 縦
@constraint(model, sum(vs[x, x] for x = 1:size) == w) # 斜め
@constraint(model, sum(vs[x, size - x + 1] for x = 1:size) == w)
# all different の制約条件
@variable(model, ns[1:size2, 1:size2], Bin)
@constraint(model, [n = 1:size2], sum(ns[n, k] for k = 1:size2) == 1)
@constraint(model, [k = 1:size2], sum(ns[n, k] for n = 1:size2) == 1)
# 盤面と数字の対応
for x = 1:size
for y = 1:size
@constraint(model, vs[x, y] == sum(ns[(x - 1) * size + y, k] * k for k = 1:size2))
end
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
for x = 1:size
for y = 1:size
print(round(Int64, value(vs[x, y])), " ")
end
println("")
end
println(value.(vs))
end
プログラムは簡単なので説明は割愛いたします。それでは実際に試してみましょう。
julia> include("magic2.jl")
solver (generic function with 1 method)
julia> @time solver(3)
結果: OPTIMAL
2 9 4
7 5 3
6 1 8
0.230510 seconds (364.44 k allocations: 17.672 MiB, 4.81% gc time, 92.96% compilation time)
julia> @time solver(4)
結果: OPTIMAL
11 5 4 14
9 7 2 16
6 10 15 3
8 12 13 1
0.095495 seconds (9.54 k allocations: 658.734 KiB)
julia> @time solver(5)
結果: OPTIMAL
14 4 13 24 10
8 22 3 9 23
21 2 12 25 5
7 19 17 6 16
15 18 20 1 11
1.864169 seconds (19.29 k allocations: 1.644 MiB)
5 行 5 列の魔方陣の場合、約 2 秒で答えを求めることができました。M.Hiroi の実行環境では、このような単純な制約では 5 行 5 列が限界だと思います。対称解のチェックを入れて解を制限すると、もう少し速くなるかもしれません。興味のある方は試してみてください。
次は下図に示す経路で、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 に行く辺を Xij, その辺のコストを Wij とすると、目的関数と制約条件は次の式で表すことができます。
スタートとゴール以外の頂点の場合、辺がつながっていないときは入力が 0 で出力が 0 になります。辺がつながっているときは入力が 1 で出力が 1 になります。これを制約条件 3 で表しています。
プログラムは隣接行列を使うと簡単です。次のリストを見てください。
リスト : 最短経路問題
using JuMP, HiGHS
#
# 辺の重み
#
# 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 = length(node) # 頂点の数
start = 1 # スタート (A)
goal = 7 # ゴール (G)
model = Model(HiGHS.Optimizer)
# 変数の生成
@variable(model, xs[1:size, 1:size], Bin)
# 目的関数
@objective(model, Min, sum(ws[i, j] * xs[i, j] for i = 1:size for j = 1:size))
# 制約条件
@constraint(model, sum(xs[start, j] for j = 1:size) == 1)
@constraint(model, sum(xs[i, goal] for i = 1:size) == 1)
for i = 1:size
if i != start && i != goal
@constraint(model, sum(xs[j, i] for j = 1:size) == sum(xs[i ,j] for j = 1:size))
end
end
# 重み 0 の辺は選択しない
for i = 1:size
for j = 1:size
if ws[i, j] == 0
@constraint(model, xs[i, j] == 0)
end
end
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
for i = 1:size
for j = 1:size
if value(xs[i, j]) == 1
println(node[i], " -> ", node[j])
end
end
end
println(objective_value(model))
辺の重み (コスト) は配列 ws に格納します。ws[i, j] が Wij を表します。0 は辺がないことを表します。変数 (Bin) は二次元配列 xs にセットします。xs[i, j] が変数 Xij を表します。目的関数と制約条件は式をそのままプログラムするだけです。あとは ws[i, j] が 0 であれば、制約条件に xs[i, j] == 0 を追加します。
それでは実行してみましょう。
結果: OPTIMAL A -> B B -> D D -> E E -> G 16.0
最短経路は A -> B -> D -> E -> G で、コストは 16 になりました。頂点の数が一番少なくなる経路は A - C - E - G ですが、これだとコストが 21 になってしまいます。ちなみに、ゴールを F にすると結果は次のようになります。
結果: OPTIMAL A -> B B -> D D -> F 6.0
A -> B -> D -> F が最短経路でコストは 6 になります。これは頂点数が一番少ない経路と一致します。