Advent Calendar 19日目の記事組合せ最適化で数コロを解く
Advent Calendar 21日目の記事組合せ最適化でひとりにしてくれを解く

これなに

ビルディングを、Pythonで組合せ最適化モデルを作って解きます。
解く楽しみは、モデル化を工夫することになります。

自分でも試してみたい人は、下記を参考にしてください。

数独を通して組合せ最適化を学ぼう

問題

各数字はそのマスに立てられるビルの高さを表します。
各横行に同じ数字は入りません。
各縦列に同じ数字は入りません。
盤面の外側の数字はその数字の書かれている場所から盤面を眺めたときに同じ横行(または縦列)に見えるビルの数を表します。

下記は、左が問題で、右が答えです。

入力パラメータ

dataにヒントが入っているとします。

python

import numpy as np
from pulp import LpProblem, lpSum, lpDot, value
from ortoolpy import addvar, addvars, addbinvars
data = """\
   33  
 ..... 
 .....5
2.....1
3.....3
 .....2
  5 13 """.splitlines()
n = len(data)-2

Pythonで解く

数理モデルを作成し、解きましょう。

変数

v：各位置がどの高さか (1)
r：各位置の高さ (2)
u：各方向別の数字のある各列ごとに (3)

制約

$v_{ij}$は、1つの高さのみ (4)
rをvで表現 (5)
縦または横に同じ高さがないこと (6)
上下左右から見たときにuの合計が「数字-1」(すなわち、高くなるときにｕ==1とします) (7)

python

m = LpProblem()
v = np.array(addbinvars(n, n, n)) # (1)
r = np.array(addvars(n, n)) # (2)
def add(m, c, p, q, y, x):
    if not c.isdigit():
        return
    k = int(c)
    u = addvars(n-1) # (3)
    m += lpSum(u) == k - 1 # (7)
    vmx = r[p,q]
    for i in range(1,n):
        vnx = r[p + y*i][q + x*i]
        m += vmx + n * u[i-1] >= vnx + 1 # (7)
        m += vmx + 1 <= vnx + n - n * u[i-1] # (7)
        vtm = addvar()
        m += vmx <= vtm # (7)
        m += vnx <= vtm # (7)
        vmx = vtm
    m += vmx <= n # (7)
for i in range(n):
    for j in range(n):
        m += lpSum(v[i,j,:]) == 1 # (4)
        m += lpDot(range(n), v[i,j]) + 1 == r[i,j] # (5)
        m += lpSum(v[i,:,j]) == 1 # (6)
        m += lpSum(v[:,i,j]) == 1 # (6)
    add(m, data[i+1][  0],   i,   0,  0,  1)
    add(m, data[i+1][n+1],   i, n-1,  0, -1)
    add(m, data[  0][i+1],   0,   i,  1,  0)
    add(m, data[n+1][i+1], n-1,   i, -1,  0)
m.solve()

結果の表示

python

print(np.vectorize(value)(r).astype(int))
>>>
[[2 5 1 3 4]
 [5 4 3 2 1]
 [4 3 2 1 5]
 [1 2 5 4 3]
 [3 1 4 5 2]]

解けていることが確認できます。

以上