ポイント
- Sympyを使った2値変数の定義とQUBOモデルの構成法
- Gurobiを用いてQUBO問題を解く方法
N-QueensパズルのQUBO化
N-QueensパズルをSymPyを使ってQUBO問題として記述し，Gurobi optimizerで解くPythonプログラムです．
$n\times n$のボードに$n$個のクイーンを配置する問題ですが，QUBO化するにあたって，$n\times n$の2値行列$X=(x_{i,j})$を用いて，$x_{i,j}=1$のときにQueenが座標$(i,j)$に置かれているとみなします．このとき$X$がN-Queensパズルとして正しい解である必要十分条件は次の通りです．
- 各行各列に1が1個ずつ(one-hot制約)，かつ
- 斜め方向の各ライン上に1が0個か1個(zero-one-hot制約)．
one-hot制約とzero-one-hot制約
$n$個の2値変数$x_0,x_1,...x_{n-1}$がone-hot制約を満たすのは，二次式
$$(x_0+x_1+\cdots+x_{n-1}-1)^2$$
が最小値0をとる場合です．zero-one-hot制約を満たすのは，二次式
$$(x_0+x_1+\cdots+x_{n-1})(x_0+x_1+\cdots+x_{n-1}-1)$$
が最小値0をとる場合です．よって，全ての行と列がone-hot制約を満たす二次式と，全ての斜めのラインがzero-one-hot制約を満たす二次式を全て足せば，N-Queensパズルを解くQUBOモデルが得られます．
Pythonプログラム
PyQUBOを使用する方法もありますが、ここでは敢えてSymPyを利用しています．2値変数を扱うためsp_Binaryクラスを定義しています．
全ての行と列がone-hotで、全ての斜めのラインがzero-one-hotであるような二次式を合計し、その結果得られた式formulaを展開し、項に分解したものがself.quboです．この式が最小値0を取る時が最適となります．定数項を除いて一次項と二次項のみとすると、定数項の符号を反転させた値が最適解self.optimalになります．行と列に関する二次式の個数は$2n$であり、この値が定数項に相当するため、self.optimalの値は$-2n$となります．この一次項と二次項をGurobi Optimizerに入力し、最適化を行っています．
#!/usr/bin/env python3
import argparse
import sympy as sp
import gurobipy as gp
class sp_Binary(sp.Symbol):
def _eval_power(self, _):
return self
class NQUEEN:
def __init__(self, n):
self.n = n
self.X = [
[sp_Binary(f"x_{i}_{j}") for j in range(self.n)] for i in range(self.n)
]
def generate(self):
sumH = [sum([self.X[i][j] for j in range(self.n)]) for i in range(self.n)]
sumV = [sum([self.X[i][j] for i in range(self.n)]) for j in range(self.n)]
A = [
[self.X[i + k][k] for k in range(self.n) if 0 <= i + k and i + k < self.n]
for i in range(-self.n, 2 * self.n)
]
B = [
[self.X[i - k][k] for k in range(self.n) if 0 <= i - k and i - k < self.n]
for i in range(-self.n, 2 * self.n)
]
sumA = [sum(x) for x in A if len(x) >= 2]
sumB = [sum(x) for x in B if len(x) >= 2]
formula = (
sum([(x - 1) ** 2 for x in sumH])
+ sum([(x - 1) ** 2 for x in sumV])
+ sum([x * (x - 1) for x in sumA])
+ sum([x * (x - 1) for x in sumB])
)
self.qubo = formula.expand().as_coefficients_dict()
def gurobi(self, timelimit):
self.model = gp.Model()
self.model.setParam("TimeLimit", timelimit)
var = {}
for i in range(self.n):
for j in range(self.n):
var[self.X[i][j]] = self.model.addVar(
vtype=gp.GRB.BINARY, name=f"x_{i}_{j}"
)
obj = []
self.optimal = -self.qubo.pop(1, 0)
for key, val in self.qubo.items():
if isinstance(key, sp.Mul):
quad = key.expand().args
obj.append(val * var[quad[0]] * var[quad[1]])
else:
obj.append(val * var[key] * var[key])
self.model.setObjective(gp.quicksum(obj), sense=gp.GRB.MINIMIZE)
self.model.optimize()
def print(self):
for i in range(self.n):
for j in range(self.n):
print(int(self.model.getVarByName(f"x_{i}_{j}").X), end="")
print()
print(f"optimal={self.optimal} obtained={int(self.model.ObjVal)}")
def main():
parser = argparse.ArgumentParser(
description="Solve N-Queens puzzle through SymPy using Gurobi"
)
parser.add_argument(
"-n", default=8, type=int, help="Value of N of N-QUEENS (default:8)"
)
parser.add_argument(
"-t", default=10, type=int, help="Time limit of Gurobi optimizer"
)
args = parser.parse_args()
nqueen = NQUEEN(args.n)
nqueen.generate()
nqueen.gurobi(args.t)
nqueen.print()
if __name__ == "__main__":
main()
実行例
Changed value of parameter TimeLimit to 10.0
Prev: inf Min: 0.0 Max: inf Default: inf
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 16 physical cores, 32 logical processors, using up to 32 threads
Optimize a model with 0 rows, 64 columns and 0 nonzeros
Model fingerprint: 0xe1261ddd
Model has 792 quadratic objective terms
Variable types: 0 continuous, 64 integer (64 binary)
Coefficient statistics:
Matrix range [0e+00, 0e+00]
Objective range [0e+00, 0e+00]
QObjective range [4e+00, 4e+00]
Bounds range [1e+00, 1e+00]
RHS range [0e+00, 0e+00]
Found heuristic solution: objective 0.0000000
Presolve time: 0.00s
Presolved: 728 rows, 792 columns, 2184 nonzeros
Variable types: 0 continuous, 792 integer (792 binary)
Root relaxation: objective -6.400000e+01, 64 iterations, 0.00 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -64.00000 0 64 0.00000 -64.00000 - - 0s
H 0 0 -10.0000000 -64.00000 540% - 0s
H 0 0 -12.0000000 -64.00000 433% - 0s
H 0 0 -14.0000000 -64.00000 357% - 0s
H 0 0 -16.0000000 -64.00000 300% - 0s
0 0 -42.66667 0 64 -16.00000 -42.66667 167% - 0s
0 0 -26.00000 0 64 -16.00000 -26.00000 62.5% - 0s
0 2 -26.00000 0 64 -16.00000 -26.00000 62.5% - 0s
Cutting planes:
Gomory: 1
MIR: 368
Zero half: 146
RLT: 289
BQP: 216
Explored 153 nodes (2925 simplex iterations) in 0.74 seconds
Thread count was 32 (of 32 available processors)
Solution count 5: -16 -14 -12 ... 0
No other solutions better than -16
Optimal solution found (tolerance 1.00e-04)
Best objective -1.600000000000e+01, best bound -1.600000000000e+01, gap 0.0000%
00010000
00000001
10000000
00100000
00000100
01000000
00000010
00001000
optimal=-16 obtained=-16
参考文献
- Shunsuke Tsukiyama, Koji Nakano, Yasuaki Ito, Takashi Yazane, Junko Yano, Takumi Kato, Shiro Ozaki, Rie Mori, Ryota Katsuki: Solving the N-Queens Puzzle by a QUBO Model with Quadratic Size. Proc. of International Symposium on Computing and Networking, 59-67, 2023. (DOI)