#用Pyomo解数独
#数独
import pyomo.environ as pe
a = [1,2,3,4,5,6,7,8,9]
c = [0,1,2]
d = [1,2,3]
m = pe.ConcreteModel()
#定义变量,一个9*9的表,每个表元素定义9个可选值
m.x = pe.Var(a,a,a,within=pe.Binary,initialize=0)
#每个元素只能有一个值
def rule_constr0(m,i,j):
return sum([m.x[i,j,k] for k in a]) == 1
m.constr0 = pe.Constraint(a,a,rule=rule_constr0)
#每列元素之和为45
def rule_constr1(m,i):
return sum([k*m.x[i,j,k] for j in a for k in a]) == 45
m.constr1 = pe.Constraint(a,rule=rule_constr1)
#每列的元素不能相同
def rule_constr11(m,i,j):
return sum([m.x[i,k,j] for k in a]) == 1
m.constr11 = pe.Constraint(a,a,rule=rule_constr11)
#每行元素之和为45
def rule_constr2(m,i):
return sum([k*m.x[j,i,k] for j in a for k in a]) == 45
m.constr2 = pe.Constraint(a,rule=rule_constr2)
#每行的元素不能相同
def rule_constr21(m,i,j):
return sum([m.x[k,i,j] for k in a]) == 1
m.constr21 = pe.Constraint(a,a,rule=rule_constr21)
#每个9格小方块数字之和为45
def rule_constr3(m,i,j):
return sum([k*m.x[3*i+t,3*j+n,k] for t in d for n in d for k in a]) == 45
m.constr3 = pe.Constraint(c,c,rule=rule_constr3)
#每个9格小方块中数字不能相同
def rule_constr31(m,i,j,k):
return sum([m.x[3*i+t,3*j+n,k] for t in d for n in d]) == 1
m.constr31 = pe.Constraint(c,c,a,rule=rule_constr31)
#目标函数,在数独中,只要变量满足约束条件就行
m.obj = pe.Objective(expr = sum([m.x[i,j,k] for i in a for j in a for k in a]),sense=pe.minimize)
opt = pe.SolverFactory('glpk')
result = opt.solve(m)
jj = []
#打印结果
for i in a:
for j in a:
for k in a:
if(pe.value(m.x[i,j,k])==1):
jj.append(k)
print(jj)
jj = []