数独をGurobiとPythonで解く

数独を0-1整数計画問題として定式化して、前回に引き続きGurobiとPythonを使って解いた。

まずは0-1整数計画問題として定式化するために数独を10進数から2進数の形式に変形する。
ある問題の解答例：
変換前

変換後

こうすると特定の行・列・ブロックにおいてそれぞれの数字が一度しか現れないことを、総和によって簡単に表現できる。
数独を解く関数sudoku_solver(...)の入出力は変換前の形式（文字列のリスト）とする。

ソースコード

itertools.productは直積を返す。for文のネストが深くならないように使ってるだけ。

#coding: UTF-8

from gurobipy import *
from itertools import product

def sudoku_solver(problem):
    m = Model("sudoku")

    # 01形式の決定変数の辞書．キーは（行，列，そのマスの数字）
    solution01 = {}
    for c,r,n in product(xrange(1,10), repeat=3):
        solution01[c,r,n] = m.addVar(vtype=GRB.BINARY,
                                     name=str(c)+str(r)+str(n))

    m.update()

    # 目的関数の設定 （実行可能解が欲しいだけなのでテキトーで）
    m.setObjective(solution01[1,1,1])
    
    # パラメータ：与えられた問題を01形式の問題に変換
    problem01 = {}
    for c,r in product(xrange(1,10), repeat=2):
        for n in xrange(1, 10):
            problem01[c,r,n] = 0
        
        if problem[c-1][r-1] != '0':
            problem01[c,r,int(problem[c-1][r-1])] = 1

    # 問題の時点で埋まっているマスはそれに決定
    for c,r,n in product(xrange(1,10), repeat=3):
        if problem01[c,r,n] == 1:
            m.addConstr(solution01[c,r,n] == 1)
                
    # すべてのマスに1から9のいずれかの数字を一つだけ割当てる
    for c,r in product(xrange(1,10), repeat=2):
        m.addConstr(sum(solution01[c,r,n] for n in xrange(1, 10)) == 1)
                
    # 各行は1から9の数字をそれぞれ1つずつ含む
    for c,n in product(xrange(1,10), repeat=2):
        m.addConstr(sum(solution01[c,r,n] for r in xrange(1, 10)) == 1)

    # 各列は1から9の数字をそれぞれ1つずつ含む
    for r,n in product(xrange(1,10), repeat=2):
        m.addConstr(sum(solution01[c,r,n] for c in xrange(1, 10)) == 1)
            
    # 各ブロックは1から9までの数字をそれぞれ1つずつ含む
    for c,r,n in product(xrange(1,10,3), xrange(1,10,3), xrange(1,10)):
        m.addConstr(sum(solution01[y,x,n] for y,x in 
                            product(xrange(c,c+3),xrange(r,r+3))) == 1)

    m.optimize()

    # 01形式の解をもとの問題の形式に変換
    solution, count = [[]], 0
    for v in m.getVars():
        if v.x == 1:
            solution[-1].append(v.varName[-1])
            count += 1
            if count%9 == 0:
                solution[-1] = ''.join(solution[-1])
                if count == 81:
                    break
                else:
                    solution.append([])

    return solution

if __name__ == '__main__':
    problem = ['003020600',
            '900305001',
            '001806400',
            '008102900',
            '700000008',
            '006708200',
            '002609500',
            '800203009',
            '005010300']
        
    for column in sudoku_solver(problem):
        print column

実行結果

Optimize a model with 356 rows, 729 columns and 2948 nonzeros
Presolve removed 356 rows and 729 columns
Presolve time: 0.00s

Explored 0 nodes (0 simplex iterations) in 0.01 seconds
Thread count was 1 (of 2 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0%
483921657
967345821
251876493
548132976
729564138
136798245
372689514
814253769
695417382

解けた。まあ、ソルバーが解いてるだけで数独を解くアルゴリズム自体は全く考えてないけど…