N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解く

# ポイント
* Sympyを使った2値変数の定義とQUBOモデルの構成法
* Gurobiを用いてQUBO問題を解く方法

# N-QueensパズルのQUBO化
[N-Queensパズル](https://ja.wikipedia.org/wiki/%E3%82%A8%E3%82%A4%E3%83%88%E3%83%BB%E3%82%AF%E3%82%A4%E3%83%BC%E3%83%B3)をSymPyを使ってQUBO問題として記述し，[Gurobi optimizer](https://www.gurobi.com/solutions/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に入力し、最適化を行っています．

```Python
#!/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](https://doi.org/10.1109/CANDAR60563.2023.00015))