数独谜题解答
本例演示如何使用混合整数线性规划（MILP）和约束编程（CP）方法求解数独谜题。
我们将使用[面向问题的优化]工作流(https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii)来解决这个问题。
安装库
如果您的环境中没有安装最新版本的JuMP 软件包，请取消注释并运行下面的对话框：
Pkg.add(["JuMP", "HiGHS"])
#Pkg.add("JuMP");
要在安装完成后启动新版本的程序库，请单击 "我的账户 "按钮：
然后点击 "停止 "按钮：
按下 "Start Engee "按钮重新启动会话：
任务描述
数独是一款逻辑组合数字排位益智游戏。游戏的目标是在 9x9 的网格中填入数字，使每列、每行和九个 3x3 的子网格中都包含 1 到 9 的所有数字。
游戏规则
1.每行必须包含从 1 到 9 的数字，不得重复。
2.每列必须包含从 1 到 9 的数字，且不重复。
3.每个 3x3 子网必须包含从 1 到 9 的数字，且不重复。
为使解题成为可能，我们提供了初始提示：
.png)
数独解题示例
.png)
数独谜题从简单到极其复杂不等，其难度通常取决于所提供的初始线索的数量以及解题策略的复杂程度。一道结构合理的数独谜题应该只有一种解法。
连接图书馆
连接图书馆JuMP ：
using JuMP;
连接求解器库HiGHS ：
using HiGHS;
连接图形库Plots ：
using Plots;
连接随机数生成库Random ：
using Random;
创建数独可视化功能
创建一个可视化数独的函数plot_sudoku() ：
function plot_sudoku(grid)
p = plot(aspect_ratio=:equal, legend=false, axis=false, ticks=false)
for i in 0:9
plot!(p, [i, i], [0, 9], color=:black, linewidth=(i % 3 == 0 ? 2 : 1))
plot!(p, [0, 9], [i, i], color=:black, linewidth=(i % 3 == 0 ? 2 : 1))
end
for i in 1:9, j in 1:9
if grid[i,j] != 0
annotate!(p, j-0.5, 9.5-i, text(string(grid[i,j]), 14))
end
end
return p
end
在本示例中，我们将使用该函数将数独可视化。
生成数独谜题
在本节中，我们将编写一个生成数独谜题的程序。
创建一个函数generate_full_grid() ，负责构建数独网格：
function generate_full_grid()
grid = zeros(Int, 9, 9)
fill_grid!(grid, 1)
return grid
end
创建一个递归函数fill_grid!() ，负责在数独网格中填入数字：
function fill_grid!(grid, pos)
if pos > 81
return true
end
row, col = divrem(pos - 1, 9) .+ 1
nums = shuffle(1:9)
for num in nums
if is_valid(grid, row, col, num)
grid[row, col] = num
if fill_grid!(grid, pos + 1)
return true
end
grid[row, col] = 0
end
end
return false
end
该函数实现了一种回溯算法，用有效数字填充数独网格。
创建一个函数is_valid() ，用于检查给定的数字 (num) 是否允许放在数独网格的某个单元格（行、列）中：
function is_valid(grid, row, col, num)
for i in 1:9
if grid[row, i] == num || grid[i, col] == num
return false
end
end
box_row, box_col = 3 * ((row - 1) ÷ 3) + 1, 3 * ((col - 1) ÷ 3) + 1
for i in 0:2, j in 0:2
if grid[box_row + i, box_col + j] == num
return false
end
end
return true
end
该函数对于在数独网格填充过程中遵循所有三条规则至关重要。
创建一个函数remove_numbers!() ，负责从填满的数独网格中删除数字，以创建一个谜题：
function remove_numbers!(grid, num_to_remove)
positions = [(i, j) for i in 1:9 for j in 1:9]
shuffle!(positions)
for _ in 1:num_to_remove
for (row, col) in positions
if grid[row, col] != 0
temp = grid[row, col]
grid[row, col] = 0
solution_count = count_solutions(grid)
if solution_count != 1
grid[row, col] = temp
else
break
end
end
end
end
end
这是制作不同难度数独谜题的关键部分，因为移除的单元格数量会影响谜题的难度。
创建一个递归函数count_solutions() ，用于计算给定数独网格的有效解数：
function count_solutions(grid, limit=2)
empty_cell = findfirst(iszero, grid)
if empty_cell === nothing
return 1
end
row, col = Tuple(empty_cell)
count = 0
for num in 1:9
if is_valid(grid, row, col, num)
grid[row, col] = num
count += count_solutions(grid, limit - count)
grid[row, col] = 0
if count >= limit
return count
end
end
end
return count
end
该函数对于确保从网格中删除数字（在函数remove_numbers! 中）不会导致多解至关重要。
创建一个函数generate_sudoku_puzzle() ，用给定的线索数生成数独：
function generate_sudoku_puzzle(num_clues)
full_grid = generate_full_grid()
num_to_remove = 81 - num_clues
remove_numbers!(full_grid, num_to_remove)
return full_grid
end
数独谜题的单解理论上最少提示数为 17。
我们编写的程序可以稳定地生成具有单一解法的数独谜题，其最小提示数介于 20 到 30 之间，这对于本例来说绰绰有余。
生成有 25 条提示的数独：
sudoku_puzzle = generate_sudoku_puzzle(25);
根据线索可视化数独
plot_sudoku(sudoku_puzzle)
混合整数线性规划问题的表述
我们将用混合整数线性规划来解决这个问题。我们的目标是将谜题建模为一个可行性问题，其中每个单元格都必须满足某些约束条件。
使用函数Model() 创建一个优化问题，并在括号中给出求解器的名称：
sudoku_milp = Model(HiGHS.Optimizer)
创建一个变量x, 表示一个三维数组, 其中i,j, 和k 的范围是 1 到 9。参数Bin 表示变量x 是二进制变量。
@variable(sudoku_milp, x[i = 1:9, j = 1:9, k = 1:9], Bin);
数组x 包括
*i （行索引）：表示数独网格中的行号，范围为 1 到 9
*j （列索引）：表示数独网格中某一列的编号，范围也是 1 到 9
*k （数字索引）：表示可放入单元格的潜在数字（从 1 到 9） (i,j)
因此，x 的二进制表示法意味着：
* - 数字 8 放在第 3 行第 3 列的单元格中
* = 0 - 数字 8 未放在该单元格中
使用循环for 循环浏览行i 和列j ，为优化问题创建以下条件：
这个条件确保数独网格中的每个单元格都包含一个介于 1 和 9 之间的数字。
for i in 1:9
for j in 1:9
@constraint(sudoku_milp, sum(x[i, j, k] for k in 1:9) == 1)
end
end
使用循环for 循环浏览行i 和列j ，为优化问题创建以下条件：
该约束确保k 的每个数字在ind 的每一行中都正好出现一次。
该约束确保k 的每个数字在ind 的每一列中恰好出现一次。
for ind in 1:9
for k in 1:9
@constraint(sudoku_milp, sum(x[ind, j, k] for j in 1:9) == 1)
@constraint(sudoku_milp, sum(x[i, ind, k] for i in 1:9) == 1)
end
end
上述单元格中的代码为模型添加了 162 个约束：
- 81 个行约束（9 行 × 9 位数）
- 81 个列约束（9 列 × 9 位数）
有了这两个条件，我们就满足了数独的前两条规则：
- 每一行必须包含从 1 到 9 的数字，且不重复
- 每列必须包含从 1 到 9 的数字，且不重复
使用循环for ，为优化问题创建以下条件：
i 的最外层循环在 1:3:7 处搜索第 1、4 和 7 行，这是每个 3x3 子网的初始行。
j 在 1:3:7 处的中间循环遍历第 1、4 和 7 列，它们是每个 3x3 子网的初始列。
k 的最内层循环在 1:9 处搜索所有可能的数字（从 1 到 9）。
通过r 和c 变量，您可以为数独中的每个 3x3 子网定义约束条件。
*r - 代表 3x3 子网格中的行索引，取值范围为i 至 。i+2
*c - 代表 3x3 子网中的列索引，取值范围为j 至j+2
for i in 1:3:7
for j in 1:3:7
for k in 1:9
@constraint(sudoku_milp, sum(x[r, c, k] for r in i:(i+2), c in j:(j+2)) == 1)
end
end
end
上述单元格中的代码为模型添加了 81 个约束条件（9 个子网格中每个网格 9 个数字），确保每个 3x3 子网格包含从 1 到 9 的所有数字，且不重复。这就是我们如何实现数独谜题的第三条规则。
利用循环for 在单元格中循环，使用函数fix() 将二进制变量x 中的值 1 (true) 放到有初始线索的单元格中 (sudoku_puzzle) ：
for i in 1:9
for j in 1:9
if sudoku_puzzle[i, j] != 0
fix(x[i, j, sudoku_puzzle[i, j]], 1; force = true)
end
end
end
您已经创建了满足所有三个解谜规则的条件，并将线索输入变量x 。现在您可以解决优化问题了：
optimize!(sudoku_milp)
输出解题结果的状态：
status = termination_status(sudoku_milp)
println("Статус: ", status)
if status == MOI.OPTIMAL || status == MOI.FEASIBLE_POINT
println("Оптимальное решение найдено")
else
println("Оптимальное решение не найдено")
end
最佳**状态意味着求解器找到了问题的全局最优解。
将二进制变量x 的值保存在变量x_val 中：
x_val = value.(x)
您可以确保解决方案中的每个值k 都已在变量x_val 中分配了网格位置。
创建一个新的 9x9 矩阵，其中填满零，用于将结果从二进制格式转换为标准数独网格格式。参数Int 表示矩阵值为整数。
sol_milp = zeros(Int, 9, 9);
使用for 循环遍历数独网格中的所有单元格和可能的数字，将三维数组x_val 的值转换为二维整数数组sol_milp 。
sol_milp = zeros(Int, 9, 9)
for i in 1:9
for j in 1:9
for k in 1:9
if round(Int, x_val[i, j, k]) == 1
sol_milp[i, j] = k
end
end
end
end
现在，变量sol_milp 包含了二维整数数组格式的实际谜题答案。
将结果可视化：
plot_sudoku(sol_milp)
我们已经演示了使用混合整数线性规划来解决问题。求解器可以找到满足数独解题规则的可能数字配置，这些规则被表述为行、列和子网格约束条件。
编程问题的约束条件表述
你也可以用约束编程来模拟这个问题。让我们应用条件AllDifferent() （"全部不同"），它规定向量中不能有两个元素取相同的值。
使用函数Model() 创建一个优化问题，并在括号中指定求解器的名称：
sudoku_cp = Model(HiGHS.Optimizer)
创建一个变量sudoku_cp ，包含一个二维数组 。每个变量的范围从 1 到 9。参数Int 表示变量的值是整数。
@variable(sudoku_cp , 1 <= x[1:9, 1:9] <= 9, Int);
设置适用于数独网格每行i 的条件。MOI.AllDifferent(9) 确保i 行中的所有元素都是不同的。
这意味着从 1 到 9 的每个数字都要在每一行中出现一次。
@constraint(sudoku_cp , [i = 1:9], x[i, :] in MOI.AllDifferent(9));
设置适用于数独网格中每一列j 的条件。MOI.AllDifferent(9) 保证j 列中的所有元素都是唯一的。
这意味着从 1 到 9 的每个数字在每一列中都正好出现一次。
@constraint(sudoku_cp, [j = 1:9], x[:, j] in MOI.AllDifferent(9));
设置适用于数独网格中每个 3x3 子网格的条件。MOI.AllDifferent(9) 确保每个子网格中的所有元素都是不同的。这意味着从 1 到 9 的每个数字都会在每个 3x3 子网格中恰好出现一次。
for i in 1:3:7, j in 1:3:7
@constraint(sudoku_cp, vec(x[i:i+2, j:j+2]) in MOI.AllDifferent(9))
end
对于[i, j] 中的每个单元格，使用函数fix 将sudoku_puzzle 中的原始线索转移到变量x 中：
for i in 1:9, j in 1:9
if sudoku_puzzle[i, j] != 0
fix(x[i, j], sudoku_puzzle[i, j]; force = true)
end
end
您已经为所有三个谜题的解题规则创建了满足条件，并将线索放入变量x 中。现在您可以解决优化问题了：
optimize!(sudoku_cp)
输出解题结果的状态：
status = termination_status(sudoku_cp)
println("Статус: ", status)
if status == MOI.OPTIMAL || status == MOI.FEASIBLE_POINT
println("Оптимальное решение найдено")
else
println("Оптимальное решение не найдено")
end
最佳**状态意味着求解器找到了问题的全局最优解。
将结果保存到变量中。得到的数值类型为Float64 。您可以使用函数round 和参数Int, 将数值转换为整数。
sol_cp = round.(Int, value.(x));
将得到的数独解法可视化：
plot_sudoku(sol_cp)
您可以比较使用混合整数线性编程和约束编程得到的解法：
sol_milp == sol_cp
两种解题方法的结果是相同的。
结论
在本例中，我们使用了两种不同的方法--混合整数线性规划（MILP）和约束规划（CP）--来解决数独谜题。这两种方法都能解决数独谜题，但它们提出和解决问题的方式却有本质区别，体现了建模和解决问题的不同特点。在解决本例中的问题时，我们可以强调这两种方法的一些不同之处。
混合整数线性规划：
** 使用二进制变量
@variable(sudoku_milp, x[i = 1:9, j = 1:9, k = 1:9], Bin)
- 使用线性约束来执行数独规则
@constraint(sudoku_milp, sum(x[i, j, k] for k at 1:9) == 1)
- 模型较大，变量较多，但约束较简单
一般对模型的微小变化更稳健，随着问题规模增大，扩展性更好
** 约束编程：**
- 使用整数变量
@variable(sudoku_cp, 1 <= x[1:9, 1:9] <= 9, Int)
- 使用全局约束，如
AllDifferent，来执行数独规则
@constraint(sudoku_cp, [i = 1:9], x[i, :] in MOI.AllDifferent(9))
- 模型更小，变量更少，但约束更复杂
可能对问题的大小和结构更敏感
在实践中，如何在混合整数线性规划和约束规划之间做出选择，往往取决于问题的具体特征，对于某些问题，一种方法比另一种方法更适合。在某些情况下，结合两种方法的混合方法可能有效。