@cleardusk 2017-07-02T15:57:55.000000Z 字数 4900 阅读 2217

# 优化角度解数独

2017寒假

## 介绍

MIP(Mixed Interger Programming) 即整数规划可以解数独，以及数独的一些变种，如对角数独、杀手数独(Killer Sudoku)

## 标准数独

### Modeling

$A \in \mathbb{Z}^{d \times d}$ 表示数独矩阵，设一个三维的变量 $x \in \mathbb{Z}^{d \times d \times d}$ (标准数独 $d=9$)，用 $x(i,j,v)$ 表示 $A(i,j)$ 是否取值 $v$，就这个地方稍微有点 trick，如果觉得有点别扭，just accept it.

$x$ 的取值为 $\{0,1\}$$1$ 表示 $(i,j)$ 元素取值 $v$$0$ 则不取。

### Solving

Mixed Integer Linear Programming problems are generally solved using a linear-programming based branch-and-bound algorithm.

from gurobipy import *import math# suppose puzzle is a nested list representing the origion sudoku problem，* will be set 0 in puzzle"""* * 5 3 * * * * * 8 * * * * * * 2 * * 7 * * 1 * 5 * * 4 * * * * 5 3 * * * 1 * * 7 * * * 6 * * 3 2 * * * 8 * * 6 * 5 * * * * 9 * * 4 * * * * 3 * * * * * * 9 7 * * """model = Model('sudoku')# 3D array, n**3 variablevars = model.addVars(n, n, n, vtype=GRB.BINARY, name='G')# Fix specified valuefor i in range(n):    for j in range(n):        if puzzle[i][j] > 0:            v = puzzle[i][j] - 1            vars[i, j, v].LB = 1  # each var in vars are BINARY, if set LB to 1, it is 1 exactly# Constraint 1: each cell must take one variablemodel.addConstrs((vars.sum(i, j, '*') == 1 for i in range(n) for j in range(n)), name='V')# Constraint 2: each value must appears once per rowmodel.addConstrs((vars.sum(i, '*', v) == 1 for i in range(n) for v in range(n)), name='R')# Constraint 3: each value must appears once per columnmodel.addConstrs((vars.sum('*', j, v) == 1 for j in range(n) for v in range(n)), name='C')# Constraint 4: each value must appears each subgrids = int(math.sqrt(n))  # subgrid sizefor v in range(n):    for i0 in range(s):        for j0 in range(s):            model.addConstr(quicksum(vars[i, j, v] for i in range(i0 * s, (i0 + 1) * s)                                     for j in range(j0 * s, (j0 + 1) * s)) == 1, name='Sub')model.optimize()# parse the result"""1 4 5 3 2 7 6 9 8 8 3 9 6 5 4 1 2 7 6 7 2 9 1 8 5 4 3 4 9 6 1 8 5 3 7 2 2 1 8 4 7 3 9 5 6 7 5 3 2 9 6 4 8 1 3 6 7 5 4 2 8 1 9 9 8 4 7 6 1 2 3 5 5 2 1 8 3 9 7 6 4"""

Gurobi 的 Python API 初步使用不太容易，语法稍微晦涩，just accept it.

### 补充

$d = 9$ 时是标准数独，$d = n^2$ 时基本都能求解的，我只试了下 $d=16$ 的，也是毫秒级别内解决。MIP 的复杂度，有兴趣可以仔细研究下。

## 对角数独

for v in range(n):    model.addConstr(quicksum(vars[i, i, v] for i in range(n)) == 1, name='D')    model.addConstr(quicksum(vars[i, n - 1 - i, v] for i in range(n)) == 1, name='D')
9 * * 1 * 5 * * 8 * * * * * * * * * * * * * 4 * * * * * * 5 * * * 8 * * * * * * * * * * * 6 9 * * * * * 3 4 * * * * * * * * * 4 7 * * 5 * * 2 3 * 6 2 7 * 1 4 8 * # 注意对角线元素9 2 4 1 3 5 6 7 8 5 1 6 2 8 7 3 4 9 7 8 3 6 4 9 1 5 2 1 3 5 4 7 2 8 9 6 2 4 8 9 6 3 5 1 7 6 9 7 5 1 8 2 3 4 8 5 9 3 2 4 7 6 1 4 7 1 8 5 6 9 2 3 3 6 2 7 9 1 4 8 5 

## 杀手数独

2 1 5 6 4 7 3 9 8 3 6 8 9 5 2 1 7 4 7 9 4 3 8 1 6 5 2 5 8 6 2 7 4 9 3 1 1 4 2 5 9 3 8 6 7 9 7 3 8 1 6 4 2 5 8 2 1 7 3 9 5 4 6 6 5 9 4 2 8 7 1 3 4 3 7 1 6 5 2 8 9 

#!/usr/bin/env python# coding: utf-8from gurobipy import *from collections import namedtuple# 4x4 demo# grid = [[1, 1, 2, 2], [1, 1, 2, 2], [3, 4, 5, 6], [3, 4, 5, 6]]# grid_sum = {1: 10, 2: 10, 3: 5, 4: 5, 5: 5, 6: 5}# 9x9 demogrid = [[ 1,  1,  2,  2,  2,  3,  4,  5,  6],        [ 7,  7,  8,  8,  3,  3,  4,  5,  6],        [ 7,  7,  9,  9,  3, 10, 11, 11,  6],        [12, 13, 13,  9, 14, 10, 11, 15,  6],        [12, 16, 16, 17, 14, 10, 15, 15, 18],        [19, 16, 20, 17, 14, 21, 22, 22, 18],        [19, 20, 20, 17, 23, 21, 21, 24, 24],        [19, 25, 26, 23, 23, 27, 27, 24, 24],        [19, 25, 26, 23, 28, 28, 28, 29, 29]]grid_sum = {1: 3, 2: 15, 3: 22, 4: 4, 5: 16, 6: 15, 7: 25, 8: 17, 9: 9, 10: 8, 11: 20, 12: 6, 13: 14, 14: 17, 15: 17,            16: 13, 17: 20, 18: 12, 19: 27, 20: 6, 21: 20, 22: 6, 23: 10, 24: 14, 25: 8, 26: 16, 27: 15, 28: 13, 29: 17}def fancy_print_sol(sol):    n = len(sol[0])    out = '\n'    for i in range(n):        for j in range(n):            out += str(sol[i][j]) + ' '        out += '\n'    print(out)def solve():    n = len(grid[0])    model = Model('sudoku_killer')    vars = model.addVars(n, n, n, vtype=GRB.BINARY, name='G')    # s.t    model.addConstrs((vars.sum(i, j, '*') == 1 for i in range(n) for j in range(n)), name='V')    model.addConstrs((vars.sum(i, '*', v) == 1 for i in range(n) for v in range(n)), name='R')    model.addConstrs((vars.sum('*', j, v) == 1 for j in range(n) for v in range(n)), name='C')    # get (i,j) pair    grid_type = {grid[i][j]: [] for i in range(n) for j in range(n)}    Pos = namedtuple('Pos', ('row', 'col'))    for i in range(n):        for j in range(n):            grid_type[grid[i][j]].append(Pos(i, j))    # s.t    for tp in grid_type.keys():        ps = grid_type.get(tp)        model.addConstr(quicksum(vars[i, j, v] * (v + 1) for v in range(n) for i, j in ps) == grid_sum.get(tp), name='S')    model.optimize()    # parse solution    solution = model.getAttr('X', vars)    sol = []    for i in range(n):        row = []        for j in range(n):            for v in range(n):                if solution[i, j, v] == 1:                    row.append(v + 1)        sol.append(row)    fancy_print_sol(sol)if __name__ == '__main__':    solve()

• 私有
• 公开
• 删除