社区首页 >专栏 >如何用模拟退火算法解数独

如何用模拟退火算法解数独

HuangWeiAI

发布于 2019-10-15 15:22:58

1.8K00

代码可运行

文章被收录于专栏：浊酒清味浊酒清味

运行总次数：0

代码可运行

数独介绍

想必大家都看过或者玩过数独游戏吧。数独游戏是源自18世纪瑞士的一种数学游戏。是一种运用纸、笔进行演算的逻辑游戏。

如上图所示，玩家需要根据9×9盘面上的已知数字，推理出所有剩余空格的数字，并满足每一行、每一列、每一个粗线宫（3*3）内的数字均含1-9，不重复：

这种游戏只需要逻辑思维能力，与数字运算无关。虽然玩法简单，但提供的数字却千变万化，所以不少教育者认为数独是锻炼脑筋的好方法。

数独解法和比赛

关于数独的解法有很多种，直观法解题技巧主要有：唯一解法、基础摒除法、区块摒除法、唯余解法、矩形摒除法、单元摒除法、余数测试法等。

基础摒除法：基础摒除法就是利用1～9的数字在每一行、每一列、每一个九宫格都只能出现一次的规则进行解题的方法。
唯一解法：当某行已填数字的宫格达到8个，那么该行剩余宫格能填的数字就只剩下那个还没出现过的数字了。成为行唯一解。
唯余解法：唯余解法就是某宫格可以添入的数已经排除了8个,那么这个宫格的数字就只能添入那个没有出现的数字。

随着数独发展，各种解法也是层出不穷，可谓是百花齐放。数独游戏也有专业的比赛，比如数独世锦赛是一种世界性的数独比赛，因为参赛选手、国家之多，是目前世界上规模最大的数独比赛。每年举办一届，比赛可谓是云集了各个国家的数独高手！比赛通过层层淘汰，最后决出冠军，冠军将会被授予“数独之王”的荣誉称号！

《最强大脑》节目也引入了数独比赛：

如何用程序解数独

但是今天，我们并不打算给大家详细介绍如何给计算机设计算法来让程序自己解数独。

我们要介绍的这个算法只需要知道数独最基本的规则：并满足每一行、每一列、每一个粗线宫内的数字均含1-9，不重复。除此之外，我们并不会人为给程序设计任何“技巧”，有种“重剑无锋，大巧不工”的感觉。有种那么如此神奇叫什么呢？它就是著名的“模拟退火(simulated annealing)”算法。

模拟退火算法是寻找一个最优解的算法。用形象的话来讲，我们有一座绵延不绝的山，最优解(global minimum)在山的最低谷：

现在有一个难题，就是我们下坡很容易，上坡很难。这样的话，我们容易陷在一个位置较高的谷底，而爬不出来，也就没办法找到最低谷了。

模拟退火可以解决上面的难题，它通过模拟物质中晶体结构的形成：物质 (例如金属) 晶格中的原子可以进入具有较低能级的状态, 或者随着温度降低而保持原位。我们从一个高温出发，这时候爬坡时一件比较容易的事情，可以避免一直限在位置较高的低谷。随着不断降温，我们越来越会走到一些低的位置，最终有一定概率可以找到最低谷。

现在还缺的是，如何把数独看成一个优化问题。其实做法比较简单。我们赋予数独一个“能量”，然后能量计算的规则便是：同一个九宫格，同一行，同一列任何两个数字如果一样那么能量就是1，如果不一样那么能量就是0。只能当总能量为0的时候，此时能量最低，而且满足数独完成条件。所以通过给与数独一个能量的概念和计算规则，我们将数独问题转换成一个寻找最低能量问题。

程序解数独

我们把上面的思路用Python实现：

import numpy as npimport randomimport timeimport matplotlib.pyplot as plt
def compute_energy(snx, pro, i, j):    en=0; en1=0    for ro in range(9):        if snx[i,j]==snx[ro,j]:            en = en  + 1        if pro==snx[ro,j]:            en1 = en1 + 1        if snx[i,j]==snx[i,ro]:            en  = en  + 1        if pro==snx[i,ro]:            en1 = en1 + 1                 for ro in range(3):        for co in range(3):            if  i<3:                x0 = 0            elif i<6:                x0 = 3            else:                x0 = 6
            if  j<3:                y0 = 0            elif j<6:                y0 = 3            else:                y0 = 6 
            if sn1[i,j]==sn1[ro+x0,co+y0]:                en = en + 1            if pro==sn1[ro+x0,co+y0]:                en1 = en1 + 1    return en, en1             
start = time.time() #------------read file----------------filename = 'exam2.txt'question = open(filename, "r")
#----------initialization--------------sn = np.zeros((9,9)) for i in range(9):    line =question.readline()    sp = line.split()    for j in range(9):        sn[i,j] = float(sp[j]) sn1 = sn.copy()sn1[sn1==0]=1
#----------------program----------------for n in range(300):    print (n)    temp = 1-n/299+0.00001    beta  = 1.0/temp#----------metroplish at T-----------       for imetro in range(200):        for i in range(9):            for j in range(9):                if sn[i,j]!=0:                    continue                en = 0; en1 = 0                pro = random.randint(1,9)                if pro==sn1[i,j]:                    continue                en, en1 = compute_energy(sn1, pro, i, j)                              if (en-3)>=en1:                    sn1[i,j] = pro                elif random.random()<np.exp((en-3-en1)*beta):                     sn1[i,j] = pro total_en = 0 for i in range(9):            for j in range(9):                en, en1 =  compute_energy(sn1, pro, i, j)                total_en = total_en + en  
if total_en-3*81 == 0:    print ('done!')    print (sn1)else:    print ('Have not found the solution.')    
end = time.time()print (str(end-start))