数理最適化・(量子)アニーリングでFGO攻略 ~ 1 v.s. 1 編 ~

「FGO、面白いですよね(廃課金感)」

# 本記事の概要

ソーシャルゲーム「Fate/Grand Order」(通称FGO)を、**最適化問題として定式化しそれを量子 or シミュレーテッドアニーリングで解いてみよう**という趣旨の記事です。

# 定式化

## バイナリ変数

以下のようにバイナリ変数を定めます。

```math
x_{k} = \{0, 1\} \\
y_{k} = \{0, 1\}
```

ここで$x$はmaster(自分)、$y$はenemy(相手)を表すバイナリです。$k$はサーヴァントの種類を表す添字で、例えばmasterがアサシンをパーティに所持している場合には$x_{殺} = 1$という風になります。

## コスト関数

### サーヴァントのクラス相性行列

FGOではサーヴァントごとにクラス相性があります。

セイバー(剣)はランサー(槍)に強い: 2.0倍のダメージ
セイバー(剣)はアーチャー(弓)に弱い: 0.5倍のダメージ

という感じです。以下にサーヴァントの攻撃/防御のときの相性表をまとめておきます。

||剣|弓|槍|騎|術|殺|狂|
|---|---|---|---|---|---|---|---|
|剣|1.0|0.5|2.0|1.0|1.0|1.0|2.0|
|弓|2.0|1.0|0.5|1.0|1.0|1.0|2.0|
|槍|0.5|2.0|1.0|1.0|1.0|1.0|2.0|
|騎|1.0|1.0|1.0|1.0|2.0|0.5|2.0|
|術|1.0|1.0|1.0|0.5|1.0|2.0|2.0|
|殺|1.0|1.0|1.0|2.0|0.5|1.0|2.0|
|狂|1.5|1.5|1.5|1.5|1.5|1.5|1.5|

行が攻撃側、列が防御側です。この表を行列$(A_{ij})$と書きましょう。

### 与えるダメージとコスト関数

masterのサーヴァントがenermyのサーヴァントに与えるダメージは以下のようになります。

```math
{\rm damage} = \sum_{i} \left\{ x_{i} \left( \sum_{j} A_{ij} y_{j} \right) \right\}
```

数理最適化は「拘束条件下においてコストをいかに最小化するか」という問題なので、コスト関数には

```math
{\rm Cost} = -{\rm damage}
```

を用いることにします。

## 制約条件

### サーヴァントは一人しか召喚できない

今回は簡単のため、masterもenemyも一人ずつしかサーヴァントを呼び出せないとします。すると上述のバイナリ変数の定義と合わせて

```math
\sum_{k} x_k = 1
```

が制約条件となります。これをQUBOの制約として入れると

```math
H_a
= \left( \sum_k x_k -1 \right)^2
```

のようになります。enemy側も同様に

```math
H_b
= \left( \sum_k y_k -1 \right)^2
```

となります。

## Total QUBO

以上を踏まえて、Total QUBOは以下のようになります。

```math
{\rm Energy} 
= -{\rm damage} + \lambda_a H_a + \lambda_b H_b
```

ここで$\lambda_a, \lambda_b$はハイパーパラメータです。

# 実際のスクリプト

## main.py

### スクリプト

```python:main.py
import display_result as dr
import make_energy as me
import make_instance as mi
import set_parameters as sp
import solve_problem as sop


if __name__ == '__main__':
    # set problem
    kind, atk_matrix = mi.make_instance()
    # set cost, constraint
    model = me.make_energy(kind=kind, atk_matrix=atk_matrix)
    # set hyper parameters
    parameters = sp.set_parameters()
    # solve with OpenJij
    solution, broken = sop.solve_problem(model=model, **parameters)
    # check broken is empty of NOT
    print(broken)
    # display simulation result
    dr.display_result(solution['x'])
    dr.display_result(solution['y'])
```

### スクリプト詳細

`make_instance`でクラス相性行列とサーヴァントの種類を指定しています。
`make_energy`でバイナリ変数$x, y$を定義し、damageと制約を計算しています。
`set_parameters`でハイパーパラメータを設定しています。
`solve_problem`でアニーリングを用いて、この問題の最適解を求めることをしています。
最後に`display_result`で自分と相手のサーヴァントを表示して、結果が正しいことを確認しています。

## make_energy.py

### スクリプト

```python:make_energy.py
import numpy as np
from pyqubo import Array, Constraint, Placeholder


def make_energy(kind, atk_matrix):
    # set placeholder (hyper parameters)
    lambda_a = Placeholder('h_a')
    lambda_b = Placeholder('h_b')
    # set binary variables
    x = Array.create('x', shape=(kind), vartype='BINARY')
    y = Array.create('y', shape=(kind), vartype='BINARY')
    # convert to numpy array
    x = np.array(x)
    y = np.array(y)
    # set cost function
    atk_cost = np.dot(x, np.dot(atk_matrix, y))
    # set one-hand constraint
    h_a = Constraint((sum(x)-1)**2, label='h_a')
    h_b = Constraint((sum(y)-1)**2, label='h_b')
    # compute total energy
    energy = -atk_cost + lambda_a * h_a + lambda_b * h_b
    # compile
    model = energy.compile()
    return model
```

### スクリプト詳細

ここでは[PyQUBO](https://pyqubo.readthedocs.io/en/latest/)を用いてQUBOの自動生成を行なっています。バイナリ変数$x, y$を生成し、それらを用いてdamageと制約を計算しています。

## solve_problem.py

### スクリプト

```python:solve_problem.py
import openjij as oj


def solve_problem(model, h_a, h_b):
    # set dictionary of hyper paramters
    feed_dict = {'h_a': h_a, 'h_b': h_b, }
    # convert to qubo
    qubo, offset = model.to_qubo(feed_dict=feed_dict)
    # solve with OpenJij (SA)
    sampler = oj.SASampler()
    response = sampler.sample_qubo(Q=qubo)
    # take mininum state
    dict_solution = dict(zip(response.indices, response.min_samples['states'][0]))
    # decode for analysis
    solution, broken, energy = model.decode_solution(dict_solution, vartype='BINARY', feed_dict=feed_dict)
    return solution, broken
```

### スクリプト詳細

ここでは[OpenJij](https://openjij.github.io/OpenJijTutorial/_build/html/ja/index.html)のSA(Simulated Annealing)を用いて、最適解の算出を行なっています。この`oj.SASampler()`の行を`SQASampler`に変えるだけでシミュレーテッド量子アニーリングでの計算を行うことも可能です。

# 結果

このスクリプトを何回か実行してみます。

```
rider
caster
```

自分がライダーで、相手がキャスターです。相手よりも有利であることがわかります。

```
archer
saber
```

自分がアーチャーで、相手がセイバーです。やはり相手よりも有利なクラスを選ぶことができています。

# 結言

このようにして、この世にあるストラテジーゲームの類は量子アニーリングでシミュレーションすることができそうな気がしてきました。
今回は1対1でしたが、これを発展させてマルチ対マルチにして、実際のFGOのシステムのシミュレーションに近づける予定です。

# 参考リンク

[組合せ最適化でゲーム理論を解く](https://qiita.com/SaitoTsutomu/items/eeb19fc086640562513e)
[ペナルティ法(罰金法)](https://qiita.com/github-nakasho/items/9003a0a4e0fb5f8d1180)