N-QueenをSA（シミュレーティド・アニーリング）で解いてみる

技術
この記事は約11分で読めます。

はじめに

　Juliaを使って、N-Queen問題をシミュレーティド・アニーリングで解いてみます。

N-Queen問題とは

　8-Queen問題というものがあります。チェス盤（8 x 8）上に8個のQueenを配置するのですが、他の駒の移動可能な位置（取られる位置）に置いてはいけないという制限があります。
　Queenは、縦横斜めの8方向に遮るものがない限り進めます。（将棋の飛車と角を合わせた動き）
　したがって、縦にもQueenは一つ、横にもQueenは一つ、斜めにもQueenは高々一つということになります。

　N-Queen問題は、これをN x Nのマスに拡張したものです。
　詳細は下記を参照してください。

SA法（シミュレーティド・アニーリング法、焼きなまし法）とは

　解空間を探索する際には、最適解に近づいたかどうかを局所的に判定しながら進める方法が一般的なのですが、SAでは必ずしも最適解に近づく変化だけを行うわけではありません。ある確率で、最適解から遠ざかる変化も許容します。これは、最適解に到達する前に局所解に陥ってしまうのを防ぐ効果があります。
　詳細は、下記を参照してください。

エネルギー関数

　解空間における状態を判定するために、最適解となる制約を考えます。その制約をエネルギー関数として計算式に書き起こします。
　エネルギー値が低いほど最適解により近い状態と考えれば理解しやすいかと思います。
　ただし、同じ問題でも、制約の記述方法は一つには定まりません。よりよいエネルギー関数を考えるのも解を早く得るために必要になります。

　今回、エネルギー関数は、次の二つの記事を参考にしました。

　最終的に次のような計算式を採用しました。最適解を求めるということは、この$H$の値が最小になる状態を求めることになります。
　この計算式（$H$）では、Queenが、縦には1個、横には1個、斜めには0個か1個、そして全体のQueenの数がN個、という条件をすべて満たした場合に、最小値0になります。

$H = H_{row} + H_{col} + H_{diag} + H_{queens}$

$H_{row} = \displaystyle \sum_{row} (\displaystyle \sum_{i∈row} x_{i}-1)^{2}$

$H_{col} = \displaystyle \sum_{col} (\displaystyle \sum_{i∈col} x_{i}-1)^{2}$

$H_{diag} = \displaystyle \sum_{diag} ((\displaystyle \sum_{i∈diag} x_{i})(\displaystyle \sum_{i∈diag} x_{i}-1))$

$H_{queens} = (\displaystyle \sum_{i∈row} \displaystyle \sum_{j∈col} q_{ij}-N)^{2}$

プログラム

# N-queen by Simulated Annealing
# シミュレーティッドアニーリング（焼きなまし法）によるN-Queen
struct Board
    size::Int64
    unit::Matrix{Int8}

    rows::Vector{Vector{Tuple{Int64,Int64}}}
    cols::Vector{Vector{Tuple{Int64,Int64}}}
    diags::Vector{Vector{Tuple{Int64,Int64}}}

    function Board(size::Int64)
        # 初期化：size個のQueenを配置しておく
        u = zeros(Int8, size, size)
        for i = 1:size
            u[i,i] = 1
        end

        # energy計算用の row, col, diagの座標一覧
        r = [[(x, y) for x = 1:size] for y = 1:size]
        c = [[(x, y) for y = 1:size] for x = 1:size]

        # rightup-diagonal
        d = [[(x, a - x + 1) for x = 1:minimum([a, size]) if a - x < size] for a = 1:(2*size-1)]

        # rightdown-diagonal
        append!(d, [[(x, x) for x = 1:size]])
        append!(d, [[(x, x - b) for x = b+1:size] for b = 1:size-1])
        append!(d, [[(y - b, y) for y = b+1:size] for b = 1:size-1])
        
        new(size, u, r, c, d)
    end
end

function energy(board::Board)
    #
    # エネルギー計算
    # 最適解の場合に、エネルギー=0となる
    #
    # row
    e = sum([(sum([board.unit[x, y] for (x, y) in row]) - 1)^2 for row in board.rows])

    # col
    e += sum([(sum([board.unit[x, y] for (x, y) in col]) - 1)^2 for col in board.cols])

    # diagonal
    e += sum([s * (s-1) for s in [sum([board.unit[x, y] for (x, y) in diag]) for diag in board.diags]])

    # queen
    e += (sum([board.unit[x, y] for x = 1:board.size for y = 1:board.size]) - board.size)^2

    return e
end

function change(board::Board, x, y)
    # 盤上のQueenのある(1)なし(0)を逆転(xor)する
    board.unit[x, y] ⊻= 1
end

function display(board::Board)
    # 表示
    println()
    for y = 1:board.size
        s = []
        for x = 1:board.size
            push!(s, board.unit[x, y])
        end
        println("[" * join(s, ", ") * "]")
    end
end

function run(board::Board, temp_start, temp_end, alpha, temp_iter)
    now_energy = energy(board)

    temperature = temp_start
    while temperature > temp_end && now_energy > 0
        for i = 1:temp_iter
            # 一か所変化させる
            rx = rand(1:board.size)
            ry = rand(1:board.size)
            change(board, rx, ry)

            # エネルギー値の計算
            new_energy = energy(board)
            # 変化前とのエネルギー値の差分
            de = new_energy - now_energy

            # eval
            if exp(- de / temperature) > rand()    # 条件 (de < 0) 含む 
                # 変化を許す
                now_energy = new_energy
                #display(borad)
                #println("energy = ", now_energy, "temp = ", temperature)

                # 解に到達したら終了
                if now_energy == 0
                    return (true, now_energy, temperature)
                end
            else
                # 変化を許さないので元に戻す
                change(board, rx, ry)
            end
        end

        # 温度を減衰
        temperature = temperature * alpha
    end

    return (false, now_energy, temperature)
end

size = 8
# 開始温度
temp_start = size*size*10
# 減衰率
alpha = 0.99
# 終了温度
temp_end = 0.02
# 温度当たりの繰り返し数
temp_iter = size * size

board = Board(size)
#display(board)

(result, ene, temp) = run(board, temp_start, temp_end, alpha, temp_iter)
display(board)
println("energy = ", ene, ", temp = ", temp)
if result
    println("OK")
end

実行例

　REPL上で実行できますし、例えば「sa.jl」ファイルに保存してコンソールから実行することもできます。
　実行すると、実行後の盤面、エネルギー値、温度が表示されます。
　8-queenなら、数秒程度で結果が出ると思います。

>julia sa.jl

[0, 0, 0, 0, 0, 0, 1, 0]
[0, 0, 1, 0, 0, 0, 0, 0]
[1, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 1, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 1, 0, 0, 0]
[0, 1, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 1, 0, 0, 0, 0]
energy = 0, temp = 0.47016434910233607
OK

　この例は、最終的にエネルギーが0になっているので、最適解が得られていることがわかります。しかし、SA法では近似解を求めているので、必ずしも最適解にならない場合も出てきます。
　その場合は、何度か繰り返すことで、「energy = 0」となるパターンを得ることができます。

備考

　エネルギー値の計算は、盤上のQueenの配置に基づいて行うのですが、計算のたびに、制約をかける座標の組み合わせを計算していました。最初、盤上の斜め方向の制約（$H_{diag}）の計算がうまくいかなかったので、デバッグのために座標位置を求める処理だけを別にしました。その仕組みをそのまま残してしまったので、盤面情報（Board）に関係ない別の情報が紛れ込んでいるように見えるかもしれません。

　また、各種パラメータは、いくつか試してそれっぽいものを設定してあります。何らかの理論的根拠に基づいているわけではありません。

　最初、SA実行中の途中経過の盤面を表示していました。ところがこの途中経過表示の時間が、SA自身の計算時間よりはるかに大きくなってしまうので、最終的な結果のみを表示することにしました。

まとめ

　上記を実現したプログラムを、下記に格納してあります。

　なお、参考までに、同じアルゴリズムを使ったPythonのプログラムも格納してあります。

参考記事

コメント

タイトルとURLをコピーしました