Gurobi学习笔记——求解数独问题
本文以Gurobi官方提供的数独案例为例，将介绍以下知识点：
- 设置变量的属性Attribute
- 如何固定变量的值
- 使用生成器添加多个约束
- quicksum() 函数的使用
设置变量的属性
Gurobi中的Var类具有多个属性(Attribute)，如LB，UB，Obj等。详情可以参见文档根目录/docs/refman.pdf
这些属性可以在创建变量时，使用关键字参数传入，也可以调用相应的方法进行相应的更改
# 创建一个目标汉中的系数为2.0，变量名称为x的0-1变量
x = m.addVar(obj = 2.0, vtype = GRB.BINARY, name = 'x')
# 在得到Var对象后，再对属性进行修改
var.setAttr(GRB.Attr.UB, 0.0)
var.setAttr("ub", 0.0)
# 可以直接用.访问其属性
var.UB = 0
固定变量的值
如果要固定某一变量的值，可以令该变量的上限和下限等于同一个常数
constant = 9
var.UB = constant
var.LB = constant
注意固定变量的值不等同于变量中Start属性。Start属性为MIP问题中的初始解，可以用一些启发式规则得到该问题的某一可行解，并以此进行计算。之后的初始解可能会发生改变。因此，应注意区分。
生成器生成多个约束
Python中的生成器(generator)对象，可以在需要的时候再对某一循环进行按需调用。详情可通过廖雪峰老师的教程了解。
在Gurobi中，生成器通常按照列生成式的方法编写，只不过将外部的方括号[]变为圆括号()
# 生成一个2*x列表 （x= 1,...,9)
>>> L = [2x for x in range(10)]
>>> L
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
>>> g = (2x for x in range(10))
>>> g
<generator object <genexpr> at 0x1022ef630>
在日常建模中，通常针对某一集合中的每个元素遍历并添加约束，如：
我们可以使用m.addConstr()，针对每一种i,j的组合，为其添加约束。
Gurobi也支持将生成器传入m.addConstrs()， 以达到批量添加的目的。
# 创建3维0-1变量
vars = m.addVars(n, n, n, vtype=GRB.BINARY, name="G")
for i in range(n):
for j in range(n):
m.addConstr(vars.sum(i,j,'*')==1)
# 等价于
m.addConstrs((
vars.sum(i,j,'*')==1
for i in range(n)
for j in range(n)))
两种写法的区别仅仅是将for循环写在外面还是使用生成器，其他方面没有区别，但前者的可读性可能会稍强些
quicksum()函数的使用
quicksum(data)是Gurobi推荐的求和函数，其执行效率高于Python内置的sum()函数。因此，在大规模添加模型时，建议优先使用quicksum()
quicksum(data)接受含有Var或者表达式(LinExpr, QuadExpr)的List对象，并将其中的所有的元素相加，生成求和表达式
expr = quicksum([2*x, 2*y+1, 4*z*z])
expr = quicksum([x, y, z])
上文中的案例，除了可以用tupledict的sum函数，也可以写作quicksum。
不过tupledict的sum方法支持通配符*，书写起来更加简便。
for i in range(n):
for j in range(n):
# 对于当前的i和j，在v维度上进行求和
m.addConstr((gp.quicksum(vars[i,j,v] for v in range(n))==1)
前面说data需要List类型的对象，此处可以理解为生成器作为参数传入list()函数中，可以一次性将生成器转化为list对象
数独案例
参考自Gurobi自带的案例文件根目录/examples/python/sudoku.py
如果要运行该案例，需要在命令行中添加数据文件根目录/examples/data/sudoku1
数独盘面是个九宫，每一宫又分为九个小格。在这八十一格中给出一定的已知数字和解题条件，利用逻辑和推理，在其他的空格上填入1-9的数字。使1-9每个数字在每一行、每一列和每一宫中都只出现一次，所以又称“九宫格”。（参考自百度百科）
本案例采用三维0-1决策变量x(i,j,v)，x(i,j,v)=1代表在第i行j列的格子上填的数字为v+1(i, j, k均从0开始计数)；否则，x(i,j,v)=0。
因此，具有以下约束：
约束1：每个格子只能有一个数
约束2： 每行元素不重复
约束3： 每列元素不重复
约束4： 每个子区域内（3*3），没有重复的元素
以下是本人增加中文注释后的代码，希望通过前面的解读，还是比较易懂的。
import gurobipy as gp
from gurobipy import GRB
import math
"""
假设数独模型是这个样子
.284763..
...839.2.
7..512.8.
..179..4.
3........
..9...1..
.5..8....
..692...5
..2645..8
"""
# 假设数据存放在同目录下的data文件夹下
f = open("./data/sudoku1")
grid = f.read().split()
n = len(grid)
s = int(math.sqrt(n))
m = gp.Model()
vars = m.addVars(n, n, n, vtype=GRB.BINARY, name="G")
# 读入数据，将已知的
for i in range(n):
for j in range(n):
# 如果该位置的数已知，则通过设置LB的方式，固定变量
if grid[i][j] != '.':
# 注意此处索引方式的不同
# grid为二维list， vars为dict
v = int(grid[i][j]) - 1
vars[i, j, v].LB = 1
# 同一个位置，只能选一个数字
for i in range(n):
for j in range(n):
m.addConstr(vars.sum(i,j,'*')==1)
# 等价于
# m.addConstrs((
# vars.sum(i,j,'*')==1
# for i in range(n)
# for j in range(n)))
# 添加行约束
# 对于每行而言，每个数字只能出现一次
for i in range(n):
for v in range(n):
m.addConstr(vars.sum(i, '*', v) == 1)
# 等价于
# m.addConstrs((
# vars.sum('*', j, v) == 1
# for i in range(n)
# for v in range(n)))
# 添加列约束
# 对于每列而言，每个数字只能出现一次
for j in range(n):
for v in range(n):
m.addConstr(vars.sum('*', j, v) == 1)
#等价于
# m.addConstrs((
# vars.sum(i, '*', v) == 1
# for j in range(n)
# for v in range(n)))
# 添加子矩阵约束
# 每个子矩阵内，数字不能重复
for i0 in range(s):
for j0 in range(s):
for v in range(n):
m.addConstr(gp.quicksum(vars[i, j, v] for i in range(
i0*s, (i0+1)*s) for j in range(j0*s, (j0+1)*s)) == 1)
# 等价于
# m.addConstrs((
# gp.quicksum(vars[i, j, v] for i in range(i0*s, (i0+1)*s) for j in range(j0*s, (j0+1)*s)) == 1
# for i0 in range(s)
# for j0 in range(s)
# for v in range(n)))
# 开始优化模型
m.optimize()
# 获取vars变量的X属性
# 获得的tupledict对象，solution
solution = m.getAttr('X', vars)
for i in range(n):
sol = ''
for j in range(n):
for v in range(n):
if solution[i, j, v] > 0.5:
sol += str(v+1)
print(sol)