0-1规划解数独

背景音乐：虚拟 - 陈粒

背景

本人一直以来比较喜欢数独，水平自认为还可以，
常用的app是Sodoku Joy，最难的Maelstrom模式也能在半小时内完成。
最近因工作需要，在看运筹优化方面的内容，
突然想到数独问题也可以用0-1规划求解，非常简单！
虽然变量和约束条件会变得很多，不过这都不是事儿。
Google一下，发现Python有个线性规划模块(Pulp)正好也提供了这方面的示例。
于是，简单地做个整理。

建模思路

1. 每个单元格有9种可能取值，每种取值要么为1（被选中）要么为0（不被选中），标准的数独是9*9的网格结构，因此构造9*9*9=729个0-1变量。
2. 目标函数：由于数独问题只要靠约束条件即可求解，因此目标函数可设为常数0。
3. 约束条件:
   1）每个单元格里，只能有一种取值被选中。
   2）每行不能有重复数字；
   3）每列不能有重复数字；
   4）每个3*3小方格不能有重复数字；

案例

拿以前新闻里说的“最难数独”为例：

Python实现

1. 首先创建一个记录了初始值的字典

# 有些单元格上已有指定的数值
init_choices = {
    ("1", "1"): "5",
    ("2", "1"): "6",
    ("4", "1"): "8",
    ("5", "1"): "4",
    ("6", "1"): "7",
    ("1", "2"): "3",
    ("3", "2"): "9",
    ("7", "2"): "6",
    ("3", "3"): "8",
    ("2", "4"): "1",
    ("5", "4"): "8",
    ("8", "4"): "4",
    ("1", "5"): "7",
    ("2", "5"): "9",
    ("4", "5"): "6",
    ("6", "5"): "2",
    ("8", "5"): "1",
    ("9", "5"): "8",
    ("2", "6"): "5",
    ("5", "6"): "3",
    ("8", "6"): "9",
    ("7", "7"): "2",
    ("3", "8"): "6",
    ("7", "8"): "8",
    ("9", "8"): "7",
    ("4", "9"): "3",
    ("5", "9"): "1",
    ("6", "9"): "6",
    ("8", "9"): "5",
}

2. 然后创建数独所需的若干常量

# 创建一个列表，包含了所有出现的数值
Sequence = ["1", "2", "3", "4", "5", "6", "7", "8", "9"]

# 单元格的可能取值
Vals = Sequence

# 行列的索引
Rows = Sequence
Cols = Sequence

# 创建一个boxes列表，记录每个3*3小方格里的，各单元格的行列索引
Boxes =[]
for i in range(3):
    for j in range(3):
        Boxes += [[(Rows[3*i+k], Cols[3*j+l]) for k in range(3) for l in range(3)]]

3. 建模

import pulp

# 创建0-1变量字典，表征指定单元格内的指定数值是否被选中
choices = pulp.LpVariable.dicts("Choice", (Vals, Rows, Cols), 0, 1, pulp.LpInteger)

# 创建模型
prob = pulp.LpProblem("Sudoku Problem", pulp.LpMinimize)

# 定义目标函数：由于数独问题只要满足所有约束条件即可求出解，所以目标函数随便指定为0
prob += 0, "Arbitrary Objective Function"

# 约束条件1：在一个单元格里，只能由一个数值被选中
for r in Rows:
    for c in Cols:
        prob += pulp.lpSum([choices[v][r][c] for v in Vals]) == 1, ""

# 约束条件2：对于一个数值，只能在一行/一列/一个box里，出现一次
for v in Vals:
    for r in Rows:
        prob += pulp.lpSum([choices[v][r][c] for c in Cols]) == 1, ""

    for c in Cols:
        prob += pulp.lpSum([choices[v][r][c] for r in Rows]) == 1, ""

    for b in Boxes:
        prob += pulp.lpSum([choices[v][r][c] for (r,c) in b]) == 1, ""

# 约束条件3：由于有些单元格上已有指定的数值，因此让它们的0-1变量取值恒为1
for (r, c), v in init_choices.items():
    prob += choices[v][r][c] == 1,""

4. 模型求解
   这里需要说明的是，如果问题出的不好，数独就可能存在多解。
   因此下面用while来循环求解：
   1）每次求得一个新的解，会加上一个新的约束条件，从而保证相同的解不会再出现。
   2）直到pulp.LpStatus[prob.status]不为Optimal，那么说明已经无解了，break。
   这里用termcolor的colored函数，来print不同颜色，区分是否为初始值。

from termcolor import colored

while True:
    # 模型循环求解
    prob.solve()
    # 打印当前的状态，并判断是否有解
    print("Status:", pulp.LpStatus[prob.status])
    if pulp.LpStatus[prob.status] == "Optimal":
        # 逐行打印结果
        for r in Rows:
            # 若当前为第1、4、7行，先打印分隔符
            if r == "1" or r == "4" or r == "7":
                print("+-------+-------+-------+")
            # 记录当前行的字符内容
            s = ''
            for c in Cols:
                for v in Vals:
                    # 如果该单元格的数值被选中，则打印
                    if pulp.value(choices[v][r][c]) == 1:
                        # 若当前为第1、4、7列，先打印分隔符
                        if c == "1" or c == "4" or c =="7":
                            s += "| "
                        # 按照是否为初始值，打印不同的颜色
                        if init_choices.get((r, c)) == v:
                            s += (colored(v, 'white', 'on_magenta') + " ")
                        else:
                            s += (v + " ")
            s += "|"
            print(s)
        print("+-------+-------+-------+")
        # 添加新的约束，保证不会返回相同的解
        prob += pulp.lpSum([choices[v][r][c] for v in Vals
                                             for r in Rows
                                             for c in Cols
                                             if pulp.value(choices[v][r][c]) == 1]) <= 80
    # 如果不能找到新的解，则结束程序
    else:
        break

5. 结果
   简书里设置颜色很麻烦，直接截图吧：

   
   
   

感兴趣的话，也可以试试注释掉一个初始解，就会返回多个可行解了。

速度测试

将上述内容封装到sudoku_problem函数中，运行

%timeit sudoku_problem(init_choices)

返回结果

655 ms ± 13.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

求解数独的时间在0.655秒左右，非常好~

后续可行的计划

1. 找个随机生成数独游戏的方法，完成“发现问题-解决问题”的闭环。
2. 简单地做个UX设计，用于交互。
3. 基于Django或Flask做个web端的界面，可以方便地在上面玩数独。
4. 加入图像识别模块，这样就可以识别截图的数独问题并求解了。