数理最適化・(量子)アニーリングでFGO攻略 ~数ターン維持する 編~

# 本記事の概要

以前紹介した

[FGO最適化 1 v.s. 1 編](https://qiita.com/github-nakasho/items/0d3cc464a6d9342f13f3)
[FGO最適化 多数 v.s. 多数 編](https://qiita.com/github-nakasho/items/07aabc5a050736196a1b)

をさらに発展させて、数ターンの攻防をシミュレートさせるために改良を行います。

# 定式化

## バイナリ変数

以下のようなバイナリ変数を用います。

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

$x$は自分のパーティを表し、$y$は敵パーティを表す変数です。$t$はターン、$m$はメンバー、$c$はサーヴァントのクラスを表す変数です。例えば自分のパーティに、1ターン目に2番目のメンバーがキャスターのサーヴァントのとき$x_{1, 2, 術}=1$という風になります。

## サーヴァントのクラスの集合

この記事から敵を倒す/敵に倒されるという**戦闘不能**の事象を表現していきましょう。そのためにサーヴァントのクラスの集合として

```math
\{戦闘不能(死), セイバー(剣), アーチャー(弓), ランサー(槍), ライダー(騎), キャスター(術), アサシン(殺), バーサーカー(狂), ルーラー(秤), アヴェンジャー(讐), ムーンキャンサー(月), アルターエゴ(分), フォーリナー(降),　シールダー(盾)\}
```

のようにしましょう。すると$x_{tm, 死}=1$は、$t$ターン目に$m$番目のメンバーが戦闘不能である、というように戦闘不能を表現することができるようになります。

## コスト関数

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

戦闘不能(死)の状態が加わったため、クラスの相性表にもこれを反映させなければなりません。それを以下にまとめました(エクストラクラスも追加してあります)。

||死|剣|弓|槍|騎|術|殺|狂|秤|讐|月|分|降|盾|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|死|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|0.0|
|剣|0.0|1.0|0.5|2.0|1.0|1.0|1.0|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|弓|0.0|2.0|1.0|0.5|1.0|1.0|1.0|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|槍|0.0|0.5|2.0|1.0|1.0|1.0|1.0|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|騎|0.0|1.0|1.0|1.0|1.0|2.0|0.5|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|術|0.0|1.0|1.0|1.0|0.5|1.0|2.0|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|殺|0.0|1.0|1.0|1.0|2.0|0.5|1.0|2.0|0.5|1.0|1.0|1.0|1.0|1.0|
|狂|0.0|1.5|1.5|1.5|1.5|1.5|1.5|1.5|1.5|1.5|1.5|1.5|0.5|1.0|
|秤|0.0|1.0|1.0|1.0|1.0|1.0|1.0|2.0|1.0|0.5|2.0|1.0|1.0|1.0|
|讐|0.0|1.0|1.0|1.0|1.0|1.0|1.0|2.0|2.0|1.0|0.5|1.0|1.0|1.0|
|月|0.0|1.0|1.0|1.0|1.0|1.0|1.0|2.0|0.5|2.0|1.0|1.0|1.0|1.0|
|分|0.0|0.5|0.5|0.5|1.5|1.5|1.5|2.0|1.0|1.0|1.0|1.0|2.0|1.0|
|降|0.0|1.0|1.0|1.0|1.0|1.0|1.0|2.0|1.0|1.0|1.0|0.5|2.0|1.0|
|盾|0.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|

行が攻撃側、列が防御側です。戦闘不能は誰にもダメージを与えることもなければ誰からもダメージを受けない空気のような存在として扱います。この表を行列$(A_{ij})$と書きましょう。

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

味方のサーヴァントが敵のパーティに、$t$ターン目に与えるダメージは以下のように計算されます。

```math
{\rm damage}_t 
= \sum_i \left( \sum_m x_{tmi} \right) \left( \sum_{j} A_{ij} \sum_{m} y_{tmj} \right)
```

よってこれを全てのターンで足し合わせた

```math
{\rm damage} 
= \sum_t {\rm damage}_t
```

が最大となるように最適化を行います。

## 制約条件

### サーヴァントは一つのクラスしか持ちえない

サーヴァントが複数のクラス(例えばセイバーとランサーのクラス)を同時に持つことはできません。よってこれを数式で表すと以下のようになります。

```math
\sum_c x_{tmc} = 1
```

これをQUBOとして書き表すと

```math
H_a = \sum_{tm} \left( \sum_c x_{tmc} -1\right)^2
```

### マスターは限られた人数のサーヴァントしか召喚できない

パーティに入れられるサーヴァントの人数を$M$とすると、この制約条件は

```math
\sum_{mc} x_{tmc} = M
```

となります。これをQUBOで書き表すと

```math
H_b 
= \sum_t \left( \sum_{mc} x_{tmc} -M \right)^2
```

### 最初に召喚されたサーヴァントの状態を維持する

1ターン目にフォーリナークラスで召喚したサーヴァントが、2ターン目に突然バーサーカーになっていてはゲームとして困ったものです。そこで任意のターンでの状態$x_{tmc}$は最初の状態$x_{0mc}$を維持するような制約を入れましょう。これを数式で表すと

```math
x_{tmc} = x_{0mc}
```

そしてこれをQUBOで制約に入れようとすると

```math
H_c 
= \sum_{tmc} (x_{tmc} - x_{0mc})^2
```

のようになります。

## Total QUBO

以上よりTotal QUBOは以下のようになります。

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

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

# 実際のスクリプト

## make_energy.py

### スクリプト

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


def make_energy(max_turn, my_member, enemy_member, kind, atk_matrix, enemy):
    # set placeholder (hyper parameters)
    lambda_a = Placeholder('h_a')
    lambda_b = Placeholder('h_b')
    lambda_c = Placeholder('h_c')
    # set binary variables
    x = Array.create('x', shape=(max_turn, my_member, kind), vartype='BINARY')
    # convert to numpy array
    x = np.array(x)
    # set cost function
    atk_cost = sum([np.dot(x[t][m], sum([np.dot(atk_matrix, enemy[e]) for e in range(enemy_member)])) for m in range(my_member) for t in range(max_turn)])
    # set my menber one-hot constraint
    h_a = Constraint(sum([(sum([x[t][m][k] for k in range(kind)])-1)**2 for t in range(max_turn) for m in range(my_member)]), label='h_a') 
    # set my member constraint
    h_b = Constraint(sum([(sum([x[t][m][k] for m in range(my_member) for k in range(kind)])-my_member)**2 for t in range(max_turn)]), label='h_b')
    # set save condition constraint
    h_c = Constraint(sum([(x[t][m][k]-x[0][m][k])**2 for t in range(max_turn) for m in range(my_member) for k in range(kind)]), label='hc')
    # compute total energy
    energy = -atk_cost + lambda_a * h_a + lambda_b * h_b + lambda_c * h_c
    # compile
    model = energy.compile()
    return model
```

### スクリプト詳細

[PyQUBO]()を用いてQUBOの自動生成を行なっています。バイナリ変数$x$を$T \times M \times K$のArrayとして生成し、それらを用いてdamageと制約を計算しています。

# 結果

## ケーススタディ: 項羽・虞美人戦

第2部3章の「人智統合真国シン」の項羽(バーサーカー)と虞美人(アサシン)に最適な3人を求め、そのパーティが3ターン維持されることを確認してみましょう。結果は以下のようになります。

```
/*****turn0*****/
-----us-----
caster
caster
caster
-----enemy-----
berserker
assassin

/*****turn1*****/
-----us-----
caster
caster
caster
-----enemy-----
berserker
assassin

/*****turn2*****/
-----us-----
caster
caster
caster
-----enemy-----
berserker
assassin
```

有利なクラスであるキャスターで3ターン戦えていることがわかります。

## ケーススタディ: ライダー・キャスター・アサシンの3すくみ

今度はライダー・キャスター・アサシンが敵として出現した場合にどのようなパーティになるかを最適化してみましょう。結果は以下のようになります。


```
/*****turn0*****/
-----us-----
berserker
berserker
berserker
-----enemy-----
rider
caster
assassin

/*****turn1*****/
-----us-----
berserker
berserker
berserker
-----enemy-----
rider
caster
assassin

/*****turn2*****/
-----us-----
berserker
berserker
berserker
-----enemy-----
rider
caster
assassin
```

**バーサーカーで3ターン殴れ**という結果になりました。

# 結言

今回はこれから色々と機能を実装する上で必要になるであろう、経過ターンの概念を盛り込みました。これをさらに発展させてターンごとにメンバーを増減させたりする予定です。