### 目次
1. はじめに
2. Spam Detector の概要
3. Spam Detector の QUBO式
4. 実行結果

## 1. はじめに

初めまして！
mktと申します。

今更ながら2021年の3月に行われた Fixstars Amplify ハッカソン（詳細はURL参照）で作ったアプリケーション「スパムコメント分類器：Spam Detector」に関する記事を投稿させていただきます。
https://amplify.fixstars.com/ja/hackathon00

Spam Detector の実装や発表スライドは以下のリポジトリに置いてます。

https://github.com/morishitamakoto0330/amplify_hackathon_sample

## 2. Spam Detector の概要


![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-2557-9989-56%2F17fdee7b8d6295?alt=media&token=3a237204-0f3a-402d-b5ac-a4f684c18302)

このアプリはスパムコメント検出をクラスタリングによって行います。
2量子ビットを用いて，one-hotベクトルとしてクラスを表現します。

図中の例では，
(0, 1)：スパムコメント
(1, 0)：通常のコメント
というクラス表現となっています。

#### Spam Detector の処理の流れ
①YouTubeSpam データセットとスパム検出を行いたいコメント群を準備
②Word Mover's Distance を計算
③QUBO式を定義
④Amplifyを用いてCMOSアニーリングマシンを実行
⑤結果を出力



## 3. Spam Detector のQUBO式

![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-2557-9989-56%2F17fdee89aa23c?alt=media&token=212ebc47-69e2-4b20-a5f3-0c03276e8f9d)

本アプリは，QUBO式中の相互作用の値に Word Mover's Distance を用いています。
この値は文章間の距離を表すもので，似た文章同士では値が小さくなります。

クラスタリングを行うように定義された QUBO式においては，相互作用の値が小さいと接続された2つの量子ビットのスピンは同じになるように収束し，相互作用の値が大きいと接続された2つの量子ビットのスピンは違うものになるように収束します。

すなわち，スパムコメント同士，通常のコメント同士でそれぞれスピンが揃うように動作します。

## 4. 実行結果

実行結果を以下に示します。


---


入力1

```
1: nice song !!!!
2: I love this music :)
```
出力1

```
スパムコメントは検出されませんでした
```


---


入力2

```
1: wGzVyTLglJnRSAP
2: How to Double Your BITCOIN in 24 Hours or Less? ANY QUESTIONS? Click Here : https://doublebtc.net
3: Great site. Chic design. I found a lot of interesting things here. Visit this beautiful site –and rate it http://xxxbukake.xyz
4: リジェクトは，辛い．
5: 一日中えぺやろうか
```

出力2

```
1: wGzVyTLglJnRSAP
4: リジェクトは，辛い．
5: 一日中えぺやろうか
```

入力1のようなポジティブなコメントについてはスパム判定されず，うまく実行ができました。
入力2でも，コメント1のような無意味な文字列をスパム判定することが出来ました。

しかし，用いたデータセット（YouTubeSpam データセット）が英語のデータセットであったため，日本語のコメント2つともがスパム判定されてしまっています。

### 今後の課題
今後の課題としては以下のようなものが挙げられます。
- 日本語に対応させる
- 入力データセットに実際のブログのコメントを使用する
- QUBO式の重みの調整を行う

以上です！
最後まで読んで頂きありがとうございます。




CMOSアニーリングマシンでスパムコメントを検出してみる

# はじめに

こんにちは！Yuki_Utaaiです。

記事をご覧の皆様の中には競技プログラミングをされたことのある方もおられるかと思います。競技プログラミングの問題の一部は組合せ最適化問題に帰着できるため、アニーリングマシン（イジングマシン）で解くことができます。ここでは競技プログラミングのいくつかの問題に対して、アニーリングマシンで求解するアプローチを行ったものをご紹介いたします。
アニーリングマシンは[Fixstars Amplify AE](https://amplify.fixstars.com/ja/)を使用しています。

競技プログラミングサイトでは直接アニーリングマシンを使用することはできないので、問題のテストケースをローカルにダウンロードし、ローカルで実行することで解けたかどうかを判断しています。


### ARC001-B

[問題はこちら](https://atcoder.jp/contests/arc001/tasks/arc001_2)

問題の概要：
　エアコンの設定温度を±10℃, ±5℃, ±1℃変更できるリモコンを使って、設定温度をA℃からB℃まで変更するために押すボタンの最小回数を求める。


普通に解く場合：
・幅優先探索やワーシャルフロイド法などを使う。
・そういった難しいアルゴリズムを知らなくても、差が10℃以下になるまでは±10℃のボタンを押して、10℃以下についてはあらかじめ埋め込んでおく。（[0,1,2,3,2,1,2,3,3,2,1]）
普通で解く場合の詳細な解答についてはAtCoderサイト等の解説ページを参照ください。


アニーリングマシンで解く場合：
以下のQUBOをマシンに投入すれば良いです。

$$
H =A\Bigl(\sum_{j=1}^{N} \sum_{i=1}^{R} d_ix_{ij}-T\Bigr)^2
+B\sum_{j=1}^{N} \Bigl(\sum_{i=1}^{R} x_{ij}-1\Bigr)^2
-C\sum_{j=1}^{N} x_{0j}
$$


定式化に関する補足

![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-9826-9189-25%2F17f35f28a5a1e6?alt=media&token=0b9cb02f-08a1-4557-8181-a2152ee27bf1)

![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-9826-9189-25%2F17f35f2d6ac318?alt=media&token=d8705cd6-0f09-4b91-b846-53917905a30f)


詳細については下記のコード及びパワーポイントをご参照ください。

[コード・パワーポイントはこちら](https://github.com/YukiUtaai/amplify-hackathon)


### ARC105-A

[問題はこちら](https://atcoder.jp/contests/arc105/tasks/arc105_a)

問題の概要：
美味しさが A,B,C,D であるような 4 枚のクッキーを持っています。そのうち1枚以上のクッキーを食べて、食べるクッキーの美味しさの総和と、残るクッキーの美味しさの総和が等しくなることがあるかどうかを判定する。

普通に解く場合：
どのクッキーを食べるかの組合せは2^4-1通りなので、選んだクッキーの美味しさの総和と残るクッキーの美味しさの総和が等しくなるかどうかを全探索すれば余裕で間に合います。
（等しくなるかどうかの判定だけを行えばいいので、ソートしてA+D=B+CもしくはA+B+C=Dになるかどうかを判定すればもっと速く解けます）

アニーリングマシンで解く場合：
以下のハミルトニアンをマシンに投入すれば良いです。
（QUBOとの違いは変数がxi∊{0,1}ではなくqi∊{-1, 1}になっているところです。今回の問題ではそうすることでハミルトニアンが非常にスッキリ記述できています。）

$$
H =(Aq_{0}+Bq_{1}+Cq_{2}+Dq_{3})^2
$$

[コードはこちら](https://github.com/YukiUtaai/amplify_arc105A)

### ABC190-C

[問題はこちら](https://atcoder.jp/contests/abc190/tasks/abc190_c)

問題の概要：
・1からNまでの番号が付いたN個のお皿がある
・M個の条件があり、i個目の条件は皿ai, 皿biの両方に（1個以上）ボールが置かれるときに満たされる
・1からKまでの番号が付いたK人の人がいて、i人目の人は皿ciか皿diのどちらか一方にボールを置く
・満たされる条件の個数は最大でいくらか


普通に解く場合：
人数Kが最大で16なので、ボールを皿ciに置くか皿diに置くか2^K通りを全探索すれば良いです。

アニーリングマシンで解く場合：
以下のQUBOをマシンに投入すれば良いです。

$$
H =-A\sum_{i=1}^{M} x_{a_i}x_{b_i}
+B\sum_{i=1}^{K} \Bigl( y_{i_0}+y_{i_1}-1\Bigr)^2
+C\sum_{i=1}^{N} \Bigl(x_{i}-\sum_{i=1}^{K}y_{皿iに対応しているもの}\Bigr)^2
$$

結構複雑になってしまっていますが、
第1項：満たされる条件の数
第2項：人がciかdiどちらか一方にのみボールを置くという制約
第3項：第1項と第2項で用いる説明変数を変えているので、それらを対応させたもの
今回の問題では、人の行動によっては皿にボールが2個以上乗る可能性があり単純な定式化が難しいためこのような操作を行っています。
第2項にてどちらか一方の皿にボールが置かれるという制約を補助変数yを用いて表現
→第3項により、それに対応する変数xの0or1を決定
→第1項により、満たされる条件の数を評価
といったようなイメージです。
第3項では皿iにボールが複数乗っているようなケースでは正の値を取るため、最小となるHが満たされる条件の数とは必ずしも一致していないことには気を付けたほうが良いです。
最小となるHをを実現した変数の組を取り出し、それを用いて満たされる条件の数を計算する必要があります。


[コードはこちら](https://github.com/YukiUtaai/amplify_abc190C)

### ABC193-F 

[問題はこちら](https://atcoder.jp/contests/abc193/tasks/abc193_f)
[コードはこちら](https://github.com/YukiUtaai/amplify_abc193F)

### ABC182-C

[問題はこちら](https://atcoder.jp/contests/abc182/tasks/abc182_c)
[コードはこちら](https://github.com/YukiUtaai/amplify_abc182C)

# 終わりに

このように、競技プログラミングの問題の中にはアニーリングマシンで解くことができるものがあります（直近ではABC240-Cも解けそうです）。チュートリアルなどでよく見られる巡回セールスマン問題やスケジューリング問題だけではなく、様々な問題に対してアニーリングマシンが活用できることを知っていただけたら幸いです。競技プログラミングの問題は想定解であれば非常に高速に解けるものなのでQUBOを立式するのは遠回りにはなってしまいますが、Nが大きくなったり問題が拡張されたり複雑化されたりした場合には従来の解法では厳密解を求めることが困難となり、アニーリングマシンが力を発揮する可能性があると思います。

競技プログラミングの問題をアニーリングマシンで解いてみよう

家庭のエネルギー利用スケジュールを最適化する

# はじめに
こんにちは！keisukeです．
今回は，QUBOモデル(特に，二次割当問題)におけるエネルギー関数を詳しく見ていきます．
「アニーリングマシンで問題を解く！」の1個前の段階のお話です．


# TL;DR
- 二次割当問題とは
- 二次割当問題をQUBOモデルに落とし込む
- 二次割当問題のエネルギー関数を展開して，**エネルギー関数の中身**を詳しく見る．


# 二次割当問題(QAP)
二次割当問題，通称QAPでは，複数の工場と，工場と同数の地区が与えられます．
任意の工場間には物流量が，任意の地区間には距離が定義されています．
二次割当問題とは，「各工場間の物流量」と「工場が割り当てられた各地区間の距離」の積の総和が最小となるような，工場の地区への割り当てを求める問題です．

ではここからは，QAPを定式化していきましょう．
$F$ 個の工場の集合を $Fac = \{f_1, f_2, ..., f_F\}$，工場を割り当てるための $F$ 個の地区の集合 $Loc = \{l_1, l_2, ..., l_F\}$ と定義します．
任意の工場間には物流量が与えられており，工場 $f_i, f_j (1 \leq i, j \leq F )$ 間の物流量を $w(f_i, f_j)$ とします．
また，任意の地区間にも距離が与えられており，工場 $f_i$ が割り当てられた地区と工場 $f_j$ が割り当てられた地区間の距離を $d(f_i, f_j)$ とします．
したがって，「各工場間の物流量」と「工場が割り当てられた各地区間の距離」の積の総和 $C$ は，次の式(1)のように表せます．
$$
C = \sum_{i=1}^{F-1} \sum_{j=1, i+1}^F w(f_i, f_j) d(f_i, f_j) \tag1
$$

なお，QAPには以下の2つの制約があります．

1. 工場配置制約：任意の工場はただ一つの地区にのみ必ず配置する
2. 地区内工場制約：任意の地区にはただ一つの工場が割り当てられる

以上がQAPの問題定義および定式化(とその制約)です．


# QUBOモデルに落とし込む
## その前に
QAPをQUBOモデルに落とし込む際，以下のようなバイナリ変数を定義します．
$$
x_{i,k} = \begin{cases}
    1 & (工場 f_i を地区 l_k に配置するとき)\\
    0 & (上記以外) \tag2
\end{cases}
$$

## コスト関数
コスト関数は，「各工場間の物流量」と「工場が配置された各地区間の距離」の積の総和です．
エネルギー関数は次のように書けます．
$$
\mathrm{H}_A = \sum_{i=1}^F \sum_{j=1}^F \sum_{k=1}^F \sum_{l=1}^F w(f_i, f_j) d(l_k, l_l) x_{i, k} x_{j, l} \tag3
$$
$\mathrm{H_A}$ の最小値は $h_a$ をとります．
$h_a$ は問題に依存します．


## 制約項1：工場配置制約
これは，$i$ 番目の工場 $f_i$ はただ一つの地区にのみ配置するという制約です．
本制約をある時刻 $i$ 番目の工場 $f_i$ について定式化すると，次の式のようになります．
$$
\sum_{k=1}^F x_{i,k} = 1~(1\leq i\leq F) \tag4
$$

全ての工場 $f_i$ が上の式(4)を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_B}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_B = \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2 \tag5
$$
$\mathrm{H_B}$ の最小値は 0 をとります．


## 制約項2：地区内工場制約
これは，$k$ 番目の地区 $l_k$ にはただ一つの工場が配置されているという制約です．
本制約を $k$ 番目の地区 $l_k$ について定式化すると，次の式のようになります．
$$
\sum_{i=1}^F x_{i,k} = 1~(1 \leq k \leq F) \tag6
$$


全ての地区 $l_k$ が上の式(6)を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_C}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_C = \sum_{k=1}^F \Biggl( \sum_{i=1}^F x_{i, k} - 1 \Biggr)^2 \tag7
$$
$\mathrm{H_C}$ の最小値は 0 をとります．


## 最終的なエネルギー関数
以上より，最終的なエネルギー関数は，正のハイパーパラメータα を用いて次のように表せます．
$$
\mathrm{H} = \mathrm{H}_A + \alpha (\mathrm{H}_B + \mathrm{H}_C) \tag8
$$
エネルギー関数 $\mathrm{H}$ は最小値 $h_a$ を取り，このとき基底解となります．
基底解が得られたときのスピンが最適解となります．
なお，エネルギー関数 $\mathrm{H}$ には定数項が出現しますが，本質的にはあってもなくても変わらないので，今回は定数項は無視します．


# エネルギー関数を展開する
再掲ですが，QAPのエネルギー関数は以下の3つから構成されています．
$
\begin{aligned}
\mathrm{H}_A &= \sum_{i=1}^F \sum_{j=1}^F \sum_{k=1}^F \sum_{l=1}^F w(f_i, f_j) d(l_k, l_l) x_{i, k} x_{j, l}\\
\mathrm{H}_B &= \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2\\
\mathrm{H}_C& = \sum_{k=1}^F \Biggl( \sum_{i=1}^F x_{i, k} - 1 \Biggr)^2
\end{aligned}
$

これらのエネルギー関数を何かしらのプログラミング言語で実装しようとしたとき，どうやって書けば良いか困ってしまうことはありませんか？
$\mathrm{H}_A$ は```for```ループを4つ使って愚直に書けば良さそうですが，$\mathrm{H}_B, \mathrm{H}_C$ はシグマの中にシグマの2乗が含まれています．これは難しそうだ…

**でも大丈夫！！！**
ここから下を読めば，一発で分かります．
ということで早速，$\mathrm{H}_B, \mathrm{H}_C$ を展開していきましょう．

## $\mathrm{H}_B$ の展開
$\mathrm{H}_B$ を展開します．
$
\begin{aligned}
\mathrm{H}_B &= \sum_{i=1}^F \Biggl( \sum_{k=1}^F x_{i, k} - 1 \Biggr)^2 \\
&= \sum_{i=1}^F \Biggl(2\sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_{i, k} x_{i, l} - \sum_{k=1}^{F} x_{i, k} + 1 \Biggr) \\
&= 2 \sum_{i=1}^F \sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_{i, k} x_{i, l} - \sum_{i=1}^{F} \sum_{k=1}^{F}x_{i, k} + F
\end{aligned}
$

なぜ1行目から2行目になるのか？
これは簡単な例で試してみると，そうなることが分かります．
ここでは，工場数 $F=4$ の場合で考えてみましょう．
(以下，2乗の部分だけを計算するので，$\mathrm{H}_B$ の一番外側のシグマ $\sum_{i=1}^{F}$ は除外します)
$
\begin{aligned}
\Biggl( \sum_{k=1}^{F=4} x_k - 1 \Biggr)^2
&= \biggl( (x_1 + x_2 + x_3 + x_4) - 1 \biggr)^2 \\
&= x_1^2 + x_2^2 + x_3^2 + x_4^2 \\
&~~~~~~~~~+2x_1x_2 + 2x_1x_3 + 2x_1x_4 \\
&~~~~~~~~~+2x_2x_3 + x_2x_4 + 2x_3x_4 \\
&~~~~~~~~~-2x_1 - 2x_2 - 2x_3 - 2x_4 + 1\\
&~~~~~(…
変数x_kはバイナリ変数なので x_k^2 = x_k となるから) \\
&= 2(x_1x_2 + x_1x_3 + x_1x_4 + x_2x_3 + x_2x_4 + x_3x_4) \\
&~~~~~~~-x_1 - x_2 - x_3 - x_4 + 1 \\
&= 2\sum_{k=1}^{4-1} \sum_{l=k+1}^{4} x_k x_l - \sum_{k=1}^{4} x_k + 1 \\
&= 2\sum_{k=1}^{F-1} \sum_{l=k+1}^{F} x_k x_l - \sum_{k=1}^{F} x_k + 1
\end{aligned}
$
ほら！これで展開ができました．
お前ホンマかいな…という人は，神ツールである
[Wolfram alpha](https://www.wolframalpha.com/input/?i=%28x_1+%2B+x_2+%2B+x_3+%2B+x_4+-+1%29%5E2&lang=ja)で計算させた結果もご覧になってください．

… というわけで，エネルギー関数に2乗が含まれている場合も，今回の展開で最終的に得られた式をプログラムに落とし込めば，うまく行きそうですね！(というかうまく行きます笑)

$\mathrm{H}_C$ の展開も，扱う文字が違うだけで内容はほぼ全く同じなので，ぜひ練習として $\mathrm{H}_C$ を展開してみてください．


# おわりに
今回は，QUBOモデルで頻繁に出てくる2乗式を詳しく見ました．
このように，実際に手を動かすことが大事となってくるので，皆さんもQUBOと格闘する時はどんどん手を動かしていきましょう！

QUBOモデルを数式からしっかり理解しよう

# はじめに
こんにちは！
普段，アニーリングマシン(厳密にはイジングマシンですが)を使って数理最適化を研究しているkeisukeと申します！
今回は初投稿となるのですが，自分が研究の過程で「これは美しい…！」と感嘆したものについて投稿したいと思います！


# TL;DR
- 巡回セールスマン問題とは？
- 巡回セールスマン問題をQUBOモデルに落とし込む
- QUBO行列はどういうものなのか，実際に覗いてみる
- QUBO行列は美しい！


# 巡回セールスマン問題(TSP)とは
巡回セールスマン問題，通称TSPとは，複数の都市が与えられており，全ての都市を一度だけ訪問して出発地点に戻ってくるような巡回路の中で，総移動距離が最小となるような都市の訪問順序を求める問題です．
TSPは組合せ最適化問題としては一番有名な問題なので，ご存知の方も多いのではないでしょうか．

ここからは，TSPを定式化していきます．
$N$ 個の都市の集合を $C'=\{c_1, c_2,...,c_N\}$ とします．
任意の都市 $c_i, c_j (1 \leq i,j \leq N)$ 間には距離が定義されており，その距離を $d(c_i,c_j)$ とします．また，$i$ 番目の都市 $c_i$，$j$ 番目の都市 $c_j$ に関するバイナリ変数を $x_{i,j}$ を以下のように定義します．
$$
x_{i,j} = \begin{cases}
    1 & (c_iの次にc_jを訪問するとき) \\
    0 & (上記以外) \tag1
\end{cases}
$$
以上より，全ての都市を一度だけ訪問して出発地点に戻って くるような巡回路の総距離 $C$ は，次の式のように表せます．
$$
C = \sum_{i=1}^{N} \sum_{j=1}^N d(c_i, c_j) x_{i, j} \tag2
$$

なお，TSPには以下の2つの制約があります．
1. 同時訪問制約：一度に訪れることのできる都市はただ一つのみ
2. 訪問回数制約：任意の都市にはただ一度しか訪問しない

以上がTSPの問題定義および定式化(とその制約)です．


# QUBOモデルに落とし込む
## その前に
今回は，QUBOがどういったモデルなのかは割愛させていただきます．
多分，僕がここで書くよりググった方が遥かに分かりやすいから…泣

TSPをQUBOモデルにマッピングする際，以下のようなバイナリ変数を定義します．
$$
x_{i,j} = \begin{cases}
    1 & (時刻 t に都市 c_i を訪問する)\\
    0 & (上記以外) \tag3
\end{cases}
$$

## コスト関数
コスト関数は，ある地点から出発して全ての都市を一度だけ訪問し，元の出発地点に戻ってくるような巡回路の総移動距離です．
エネルギー関数は次のように書けます．
$$
\mathrm{H}_A = \sum_{t=1}^N \sum_{i=1}^N \sum_{j=1}^N d(c_i, c_j) x_{t, i} x_{t+1, j} \tag4
$$
$\mathrm{H_A}$ の最小値は $h_a$ をとります．
$h_a$ は問題に依存します．


## 制約項1：同時訪問制約
これは，時刻 $t$ に訪問する都市はただ一つのみであるという制約です．
本制約をある時刻 $t$ について定式化すると，次の式のようになります．
$$
\sum_{i=1}^N x_{t,i} = 1~(1\leq t\leq N) \tag5
$$

全ての時刻 $t$ が上の式を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_B}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_B = \sum_{t=1}^N ( \sum_{i=1}^N x_{t, i} - 1)^2 \tag6
$$
$\mathrm{H_B}$ の最小値は 0 をとります．


## 制約項2：訪問回数制約
これは，$i$ 番目の都市 $c_i$ はただ一度しか訪問しないという制約です．
本制約を $i$ 番目の都市 $c_i$ について定式化すると，次の式のようになります．
$$
\sum_{t=1}^N x_{t,i} = 1~(1 \leq i \leq N) \tag7
$$

全ての都市 $c_i$ が上の式を満たすときに，最小値をとるようなエネルギー関数 $\mathrm{H_C}$ を導入すると，エネルギー関数は次のように書けます．
$$
\mathrm{H}_C = \sum_{i=1}^N ( \sum_{t=1}^N x_{t, i} - 1)^2 \tag8
$$
$\mathrm{H_C}$ の最小値は 0 をとります．


## 最終的なエネルギー関数
以上より，最終的なエネルギー関数は，正のハイパーパラメータα を用いて次のように表せます．
$$
\mathrm{H} = \mathrm{H}_A + \alpha (\mathrm{H}_B + \mathrm{H}_C) \tag9
$$
エネルギー関数 $\mathrm{H}$ は最小値 $h_a$ を取り，このとき基底解となります．
基底解が得られたときのスピンが最適解となります．
なお，エネルギー関数 $\mathrm{H}$ には定数項が出現しますが，本質的にはあってもなくても変わりません．


# QUBO行列の中身を見ていく
さて，最終的なエネルギー関数 $\mathrm{H}$ をQUBO行列として取得します．
QUBO行列の中身はどのようになっているのか見ていきましょう．
今回は，都市数 $N=5$，全都市間の距離 $d(c_i, c_j)=3$，ハイパーパラメータ $\alpha=1$ の場合でやってみます．

```
import numpy as np
import plotly.express as px

class TSP:
    def __init__(self):
        self.NUM_CITY = 5 # 都市数
        self.NUM_SPIN = self.NUM_CITY * self.NUM_CITY # QUBO行列に必要なスピン数
        self.distance_matrix = np.array([[0,3,3,3,3],[3,0,3,3,3],[3,3,0,3,3],[3,3,3,0,3],[3,3,3,3,0]]) # 都市間の距離行列
        self.ALPHA = 1 # ハイパーパラメータ α

        # 各スピンに番号を割り振る
        self.spin_index = np.arange(self.NUM_CITY * self.NUM_CITY).reshape(self.NUM_CITY, self.NUM_CITY)

        # コスト関数用のQUBO行列
        self.cost_weight = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        self.cost_bias = np.zeros((self.NUM_SPIN))
        self.cost_const = np.zeros((1))

        # 制約項のQUBO行列
        self.penalty_weight = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        self.penalty_bias = np.zeros((self.NUM_SPIN))
        self.penalty_const = np.zeros((1))

    # QUBO行列を作る
    def make_qubo_matrix(self):
        # ========== ここから式(4) ==========
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY):
                for v in range(self.NUM_CITY):
                    if t < self.NUM_CITY-1:
                        idx_i = self.spin_index[t, u]
                        idx_j = self.spin_index[t+1, v]
                    elif t == self.NUM_CITY-1:
                        idx_i = self.spin_index[t, u]
                        idx_j = self.spin_index[0, v]
                    coef  = self.distance_matrix[u, v]
                    if coef == 0:
                        continue
                    self.cost_weight[idx_i, idx_j] += coef
                    self.cost_weight[idx_j, idx_i] = self.cost_weight[idx_i, idx_j]

        # ========== ここから式(6) ==========
        # 二次の項
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY-1):
                for v in range(u+1, self.NUM_CITY):
                    idx_i = self.spin_index[t, u]
                    idx_j = self.spin_index[t, v]
                    coef = 2
                    self.penalty_weight[idx_i, idx_j] += self.ALPHA * coef
                    self.penalty_weight[idx_j, idx_i] = self.penalty_weight[idx_i, idx_j]

        # 一次の項
        for t in range(self.NUM_CITY):
            for u in range(self.NUM_CITY):
                idx = self.spin_index[t, u]
                coef = -1
                self.penalty_bias[idx] += self.ALPHA * coef

        # 定数項
        self.penalty_const[0] += self.ALPHA * self.NUM_CITY

        # ========== ここから式(8) ==========
        # 二次の項
        for u in range(self.NUM_CITY):
            for t in range(self.NUM_CITY-1):
                for tt in range(t+1, self.NUM_CITY):
                    idx_i = self.spin_index[t, u]
                    idx_j = self.spin_index[tt, u]
                    coef = 2
                    self.penalty_weight[idx_i, idx_j] += self.ALPHA * coef
                    self.penalty_weight[idx_j, idx_i] = self.penalty_weight[idx_i, idx_j]

        # 一次の項
        for u in range(self.NUM_CITY):
            for t in range(self.NUM_CITY):
                idx = self.spin_index[t, u]
                coef = -1
                self.penalty_bias[idx] += self.ALPHA * coef

        # 定数項
        self.penalty_const[0] += self.ALPHA * self.NUM_CITY

        # 式(9)のHのQUBO行列を作る
        qubo_matrix = np.zeros((self.NUM_SPIN, self.NUM_SPIN))
        qubo_matrix = self.cost_weight + self.penalty_weight
        for i in range(self.NUM_SPIN):
            qubo_matrix[i, i] = self.cost_bias[i] + self.penalty_bias[i]
        fig = px.imshow(qubo_matrix)
        fig.show()
```

上記のコードを実行すると，下記のような結果が得られます．
(ちなみにplotlyというpythonライブラリ，ものすごく便利なのでぜひご活用ください！)
各セル(マス？)は各スピン間の相互作用に対応しています．
今回は都市数が5ということで，QUBOモデルに必要なスピン数は $5 \times 5=25$ となります．
ですので，QUBO行列は $25 \times 25$ の行列となるわけです．
なお，下図では色と各値が対応しており，その関係は以下のようになっています．
- 紺色：-2
- 紫色：0
- 黄色：2
- オレンジ色：3

![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-0061-7526-97%2F17edf766f022ff?alt=media&token=090f003c-eabf-4de1-aa01-c10887d820e0)

QUBO行列は対称行列となっていることが分かりますね．
これは当然といえば当然です．
なぜなら，スピン $i$ から見たスピン $j$ との間の相互作用と，スピン $j$ から見たスピン $i$ との間の相互作用は同じになるからです(逆に違ったらおかしいですよね)．
それにしてもQUBO，う，美しい…！！！


# おわりに
いかがでしたでしょうか？
QUBOは対称行列となっており，それを可視化すると美しいものになるということについて投稿させていただきました．
皆さんも，TSPに限らずいろんなQUBOモデルを可視化してみて，その美しさを味わってみてください！

QUBO行列の中身を見てみよう

巡回セールスマン問題(TSP)とは

# 概要
こんにちは、空と申します。

1年ほど前から某個別指導塾でアルバイトをしているのですが、働いていてずっと不思議に思っていたことがありました。

なんと、各コマの各生徒→講師の割り当て(以下**配置**)を手動で行っていたのです！

講師の能力、生徒の特性など様々な考慮すべき条件はあるものの、手動で配置を組むのはとても大変なので、QUBOを用いて定式化することで自動化できないかと考えました。

本記事では、この定式化と、Amplifyでの実装について記述します。

# 配置効率について
私の塾では、1人の講師に対して生徒を2人まで配置することができます。また、なるべくこの講師1生徒2の配置(以下1:2配置)を増やすように命じられています。理由としては以下の2点が挙げられます。
1. 講師1生徒1(以下1:1配置)の配置を多くしすぎると出勤講師数が多くなりより人件費がかかる。
e.g.) 生徒4人出席するとする
全て1:2配置→出勤講師数2人
全て1:1配置→出勤講師数**4人**
2. 1:1配置での授業を希望する生徒の存在。
中には、「1:1配置での授業をしてほしい」 という希望を持つ生徒もいます。このような生徒の希望に対しては、月謝を高くする形で応えています。通常生徒の1:1配置を増やしすぎてしまうと、高い月謝を払っている1:1配置希望生徒に対して不公平になってしまいます。

上記の理由から、1:2配置と1:1配置の比率がどのようになっているのかを考慮する必要があります。これを評価する数値指標として**配置効率**を定義します。配置効率は以下のような式で表されます。
$$
配置効率=\frac{実際の講師数}{必要講師数} \times {100}
$$
今回は簡単のため、1:1配置希望生徒は考慮しないこととします。そのため、必要講師数は生徒数の半分です。全て1:2配置だと配置効率は100%となり、これが最も良い値です。逆に全て1:1配置だと200%となり、これが最も悪い値となります。今回はこの配置効率を目的関数として**最小化**することを目指します。

# 制約条件
上記の配置効率に加え、今回は以下の3つの制約を組み込みます。

- 生徒制約
- 講師制約
- 生徒講師共通制約

生徒制約は、「生徒1人に対して講師はちょうど1人」という個別指導塾として当たり前の制約です。講師制約は、「1人の講師に配置できる生徒数は0人以上$F$人以下」という制約です。生徒講師共通制約は、例えば数学しか教えることのできない先生が国語の授業を担当してはならないという制約で、これも当たり前の制約です。

# QUBO変数の設定
QUBO変数を$q_{i,j}$とします。$i$が講師、$j$が生徒を表します。スピン数は$N_{\mathrm{T}}$を講師数、$N_{\mathrm{S}}$を生徒数としたとき、$N_{\mathrm{T}} \times \left(N_{\mathrm{S}}+F+1\right)$となります。QUBO変数の設定は以下です。
- $q_{i,j}=1$ $\Rightarrow$ 講師 $i$ が生徒 $j$ の授業を担当する。
- $q_{i,j}=0$ $\Rightarrow$ 講師 $i$ が生徒 $j$ の授業を担当しない。

# QUBOでの定式化
## 事前に与えるデータ
- 講師データ

講師データは講師数$\times$教科数の行列で、各講師がそれぞれどの科目の授業の担当が可能かを表しています。例えば、以下のような行列を考えます。

| 講師\教科 | 英語 | 数学 | 国語 |
| --- | --- | --- | --- |
| 0 | 1 | 0 | 0 |
| 1 | 1 | 1 | 0 |
| 2 | 0 | 1 | 0 |
| 3 | 0 | 0 | 1 |

例えば講師2の行を見ると、英語と数学の列が1となっているので、この講師2は英語と数学の授業が可能である講師だとわかります。
- 生徒データ

生徒データは生徒数$\times$教科数の行列で、各生徒が何の授業を受けるのかを表しています。例えば、以下のような行列を考えます。

| 生徒\教科 | 英語 | 数学 | 国語 |
| --- | --- | --- | --- |
| 0 | 1 | 0 | 0 |
| 1 | 0 | 1 | 0 |
| 2 | 0 | 0 | 1 |
| 3 | 1 | 0 | 0 |

例えば生徒2の行を見ると、数学の列が1となっているので、この生徒2は数学の授業を受ける生徒だとわかります。

## 制約
制約条件の項で述べた3つの制約を以下の通り定式化します。

- 生徒制約
$$
H_\mathrm{A}=
\sum_{j=0}^{N_{\mathrm{S}}-1}{\left(\sum_{i=0}^{N_{\mathrm{T}}-1}{q_{i,j}-1}\right)^{2}}
$$
- 講師制約
$$
H_\mathrm{B} =
\sum_{i=0}^{N_{\mathrm{S}}-1}
{\left(\sum_{j=0}^{N_\mathrm{T}-1}{q_{i,j}}-\sum_{k=0}^{F-1}{y_{i,k}}\right)^{2}} \\
y_{i,k}:補助変数\\
$$
実装の際はamplifyの`less_equal`を利用しました。
- 生徒講師共通制約
$$
H_\mathrm{C}=\sum_{\left(i,j\right) \in X}{q_{i,j}}\\
X:授業を担当できない(講師,生徒)のペアのリスト
$$

生徒講師共通制約のリスト$X$は、講師データと生徒データの転置の行列積を計算し、結果が0であった行と列のペアを格納することで構成します。

## 目的関数
目的関数を配置効率としたかったのですが、配置効率の定義式をそのままamplifyに落とすのは難しく、実現できませんでした。そこで、代替案として**各行のスピン和の分散**を考えました。例を示しておきます。スピン行列が以下のようであったとします。

| 講師\生徒 | 0| 1 | 2 | 3 | 4 | 5 | 6 | 7 | SUM |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 2 |
| 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 |
| 2 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 2 |
| 3 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| 4 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 2 |

この行列のSUMの列の平均を計算します。$i$行のスピン和を$Sum_i$とすると、平均$\mu$は以下のようになります。
$$
\mu = \frac{1}{N_{\mathrm{T}}}\sum_{i=0}^{N_{\mathrm{T}}-1}{Sum_i}
=\frac{1}{N_{\mathrm{T}}}\sum_{i=0}^{N_{\mathrm{T}}-1}{\sum_{j=0}^{N_{\mathrm{S}}-1}{q_{i,j}}}
$$
上記の例であれば$\mu=1.6$となります。この$\mu$を用いてSUM列のデータの分散を以下のように計算します。分散は
$$
s^2=\overline{x^2}-\left({\overline x}\right)^2
$$
で計算できるので、この式を用いると分散$H_\mathrm{F}$は、
$$
H_\mathrm{F}
=\frac{1}{N_{\mathrm{T}}}
\sum_{i=0}^{N_{\mathrm{T}}-1}
{\left(\sum_{j=0}^{N_{\mathrm{S}}-1}{q_{i,j}}\right)^2}-\mu^2
$$
となります。上記の例であれば$H_\mathrm{F}=2.8-2.56=0.24$となります。$H_\mathrm{F}$はSUM列のデータに0と2が多いほど大きくなり、1が多いほど小さくなるので、1:2配置が多くなればなるほど値が大きくなります。したがって、$H_\mathrm{F}$を**最大化**する解が最適解であるということができます。今回はこれを目的関数としました。

## エネルギー関数
- 目的関数(最大化問題$\rightarrow$負にして最小化問題に)
$$
f=-H_\mathrm{F}
$$
- 制約
$$
g=H_\mathrm{A}+H_\mathrm{B}+H_\mathrm{C}
$$

エネルギー関数は以下の通りです。今回は特にハイパーパラメータの設定は行いませんでした。
$$
H=f+g
$$

# 実装
上記の定式化を、以下のようにPythonで実装しました。
```Python
import amplify
from amplify import BinaryPoly, SymbolGenerator, sum_poly
from amplify.constraint import equal_to, less_equal

gen = SymbolGenerator(BinaryPoly)
L = gen.array(num_teacher, num_student) # QUBO変数

# 制約
H_A = sum([equal_to(sum_poly(num_teacher, lambda i: L[i, j]), 1) for j in range(num_student)]) # 生徒制約
H_B = sum([less_equal(sum_poly(num_student, lambda j: L[i, j]), 2) for i in range(num_teacher)]) # 講師制約
H_C = sum([equal_to(sum_poly(len(not_list), lambda i: L[not_list[i][0], not_list[i][1]]), 0)]) # 教務可能科目を超える授業不可

# 目的関数
inv_num_teacher = 1 / num_teacher
mu = inv_num_teacher * (sum_poly(num_teacher, lambda i: (sum_poly(num_student, lambda j: L[i, j])))) # 各行の和の平均
H_F = inv_num_teacher * (sum_poly(num_teacher, lambda i: (sum_poly(num_student, lambda j: L[i, j]) ** 2))) - mu ** 2 # 各行の和の分散

#エネルギー関数
f = -H_F
g = H_A + H_B + H_C
H = f + g
```

# おわりに
読んでいただきありがとうございました！
今回作ったものは教室業務の基本部分にしか目を向けていないので、1:1配置を希望する生徒を1:1配置するような制約や登録済みの配置の反映なども入れられたらいいと思います。githubに全体のコードを載せたのでもし興味があればご覧ください！
https://github.com/sora-tt/placement_optimizing_of_cram_school

個別指導塾の生徒配置をQUBOで自動化する

stdout

Rectangle Packing Problem (矩形パッキング問題)

## TL;DR

Amplifyはイジングマシンを手軽に扱うことのできるSDKです。

イジングマシンを使用するとき、古典マシン(一般的なCPU)でのシミュレーションと性能を比較したいケースが多々あります。
構築したイジングモデルに対してシミュレーテッドアニーリング (SA)でさくっと最適化を試したいなどが想定されます。

D-Wave Ocean SDKやPyQUBO（内部的にはdwave-neal)、OpenJijには、イジングモデル/QUBOに対するSAのソルバーが実装されています。

しかし、Amplifyでは特にSAなどの古典ソルバーは実装されていません。
（無料でイジングマシンが使用できるのでなくても全く問題ないですが）

そこで、今回はAmplifyで構築したモデルをOpenJijの古典ソルバーに投入する流れをまとめます。

## Amplifyのモデル

Amplifyでは多項式や論理式の数式オブジェクトを用いて記述し、入力模型を構築していきます。

```python
from amplify import BinaryPoly, SymbolGenerator

gen = SymbolGenerator(BinaryPoly) 
q = gen.array(2)
f = 1 - q[0] * q[1]
```

その後、構築した入力模型をイジングマシンが取り扱うことが可能な「論理模型」として抽象化します。

```python
from amplify import BinaryQuadraticModel

model = BinaryQuadraticModel(f)
```

論理模型への変換には`BinaryQuadraticModel` `IsingQuadraticModel`といった[クラス](https://amplify.fixstars.com/ja/docs/tutorial.html#id7)を使用します。

イジングマシンを利用する場合は、この論理模型をイジングマシンに投入することになりますが、今回はOpenJijのソルバーに投入するために変換していきます。


## OpenJij用の変換

OpenJijのインターフェースはイジングモデルの相互作用 $J$ と縦磁場 $h$ のリスト（QUBOの場合はQUBO行列）を渡すものになっています。

以下の関数でAmplifyの論理模型から必要な形式に変換します。

```python
from amplify import BinaryQuadraticModel, IsingQuadraticModel

def amplify_bqm_to_qubo_dict(model: IsingQuadraticModel):
    Q = {}
    const = 0
    for k, v in model.logical_model_poly.asdict().items():
        if len(k) == 1:
            Q[(k[0], k[0])] = v
        elif len(k) == 2:
            Q[k] = v
        else:
            const = v
    return Q, const

def amplify_iqm_to_ising_dict(model: BinaryQuadraticModel):
    h = {}
    J = {}
    const = 0
    for k, v in model.logical_model_poly.asdict().items():
        if len(k) == 1:
            h[k[0]] = v
        elif len(k) == 2:
            J[k] = v
        else:
            const = v
    return h, J, const
```

実際に変換関数を用いてOpenJijでSAを実行する流れは以下のようになります。

```
import openjij as oj
from dataclasses import dataclass

@dataclass
class Solution():
    # amplify.SolverSolution と同様のインターフェース
    energy: float
    frequency: int
    is_feasible: bool
    values: Dict[int, int]

sampler = oj.SASampler()
params = {'num_reads': 10}

# Ising/QUBOを変換してSASamplerで計算実行
if isinstance(model, BinaryQuadraticModel):
    Q, const = self.amplify_bqm_to_qubo_dict(model)
    result = sampler.sample_qubo(Q, **params)
elif isinstance(model, IsingQuadraticModel):
    h, J, const = self.amplify_iqm_to_ising_dict(model)
    result = sampler.sample_ising(h, J, **params)

# 結果の解析
solutions = []
for k, s in enumerate(result.states):
    # 制約違反をチェック
    c_status = list(filter(lambda c: not c[1], model.check_constraints(s)))
    values = {}
    for i in result.indices:
        values[i] = s[i]
    solution = Solution(
        energy=result.energies[k] + const,
        frequency=1,
        is_feasible=len(c_status) == 0,
        values=values,
    )
    solutions.append(solution)

print(solutions)
```

ここではAmplifyでイジングマシンを利用した場合に得られる結果のクラス `SolverSolution` と同等のインターフェースになるように `Solution` クラスを定義して利用しています。

以上でAmplifyで構築したモデルに対して手軽にSAでシミュレーションすることができます。




AmplifyのモデルをOpenJijのSASamplerに投入する

# はじめに
初めましてtyと申します。「学校のクラス分け」って聞くと、どきどきした甘酸っぱい幼少？青春？時代を思い出す人もいるのではないでしょうか。

そのクラス分けを行うのは教員なわけですが、実は相当難しいことを行っているのです。

例えば100人を3クラスに分けるとき、その組合せの総数はなんと$10^{26}$もあります。それを毎年行なっている教員ってすごい！

日々大変な教員方の負担を少しでも削減したいと思い、私たちはクラス分けを自動化するアプリケーションを製作しました。本記事ではその定式化およびAmplifyでの実装に絞って紹介します！ 

※ある程度イジングマシンや量子アニーリングをわかっている方向けの記事となっております。もしQUBOやイジングモデルと聞いてピンと来ない方は一度他のサイトや書籍で学ばれることをお勧めいたしますm(_ _)m

# クラス分けに考慮される事項
クラス分けで満足させたい事項は色々あると思います。例えば、、
- 男女比を一緒にしたい
- クラス間の学力を同じにしたい
- クラス間の居住地域の偏りを少なくしたい
- 仲の悪い生徒は別のクラスにしたい
- 運動能力の差をなくしたい
- 学力の差をなくしたい
- 特別な支援を要する生徒を振り分けたい
- ピアノを弾ける生徒を振り分けたい
- 相性の悪い親は別のクラスにしたい

ちなみに最後のは実際にあるそうm(_ _)m

# QUBO変数の設定
上記の事項をQUBOに落とし込む分けですが、QUBO変数は {0, 1} のどちらかを取る変数です。2つの状態しか表せないためクラス数が3以上の場合を考慮し工夫する必要があります。

今回は生徒一人につきクラス数分のQUBO変数を用意します。このQUBO変数を$n_{i,a}$と表記するものとします。$i$ は生徒、$a$ はクラスを表します。

- $n_{i,a}=1$ のとき生徒 $i$ がクラス $a$ に属する
- $n_{i,a}=0$ のとき生徒 $i$ がクラス $a$ に属さない

とすることでクラス数が3以上の場合も対応できます。

# QUBOの定式化
今回は以下の6つの制約を導入することにします。
  - クラス重複禁止制約
  - クラス人数均等化制約
  - 男女均等化制約
  - 学力均等化制約
  - 地域分散化制約
  - 不仲最小化制約

クラス重複禁止制約は一人が複数のクラスに属したり、どこのクラスにも属さなかったりするのを避けるための制約です。

早速ですがこれらを定式化すると以下のようになります。
$s$は生徒数、$c$はクラス数です。
![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-2915-2635-60%2F17ec956542b6d?alt=media&token=dd1081c0-8b3d-4336-bcbb-fbcbc444bcf8)

なお右側の定式化に使用する値は予め元データから抽出・整形しておきます。

# Amplifyにおけるコーディング
上記の定式化を実装します．使用言語はPythonです。

Amplifyライブラリを使用します。Amplifyライブラリの使い方は公式サイトにありますのでそちらをご覧ください。

```python
from amplify import BinaryPoly, gen_symbols, sum_poly, pair_sum

q = gen_symbols(BinaryPoly, class_num, student_num)

HA = sum([equal_to(sum_poly([q[i][j] for i in range(class_num)]), 1) for j in range(student_num)])
HB = sum_poly(class_num, lambda i: (sum_poly(student_num, lambda j: q[i][j]) - student_num_average) ** 2)
HC = sum_poly(class_num, lambda i: (sum_poly(student_num, lambda j: sex[j] * q[i][j]) ** 2))
HD = sum_poly(class_num, lambda i: (sum_poly(student_num, lambda j: grades[j] * q[i][j]) - grades_average) ** 2)
HE = sum_poly(class_num, lambda i: (pair_sum(student_num, lambda j, k: area[j][k] * q[i][j] * q[i][k])))
HF = sum_poly(class_num, lambda i: (pair_sum(student_num, lambda j, k: dislike[j][k] * q[i][j] * q[i][k])))

# ハイパーパラメータ
aA, aB, aC, aD, aE, AF = 20, 1, 1, 10, 5, 20
# H = エネルギー関数
H = aB*HB + aC*HC + aD*HD + aE*HE + aF*HF   
H = H + aA*HA

# ※amplifyの仕様上 H=HA+HB+HC+HD+HE+HF はできない模様．HAはBinaryConstraintsクラスで，HB~HFはBinaryPolyクラス．
```

こんな感じでかけました。シグマが入り乱れた制約を一行で書けるのすごいです！

なお各制約のハイパーパラメータは試行錯誤して仮決めしました。この値によってどの制約を重視するか調節できます。

# おわりに
簡単にですが定式化とAmplifyを用いた実装を紹介しました。イジングマシンや量子アニーリングの実社会への応用が進むといいですね！


学校のクラス分け自動化におけるQUBO定式化

# はじめに
Everybodyこんにちは． 深水と言います． 
この記事では，みなさんご存知の大関先生や東北大学が主催となって開催した，量子アニーリングソリューションコンテストで発表した内容について紹介させてもらいます．下のwebサイトに行けば，コンテストの概要や，全ての発表動画を見ることができます． (ちなみにまだ結果は出てませんが，優勝者はカナダのD-Wave本社まで連れて行ってくれるそうです．ファンタスティック!!)


http://www.tfc.tohoku.ac.jp/special/qca/20211218.html


私たちは，チーム名：**QuantumPredator** ，発表タイトル： **EVerify** として発表しました．
内容は **量子アニーリングを用いたEVシェアの最適化** についてです．

EVerifyの発表動画リンク↓

https://www.youtube.com/watch?v=6er7f41l-0g

チームメンバーは東北大学情報科学研究科 小林・佐藤研究室の以下３名で，この記事も共同執筆です．
東北大学情報科学研究科　博士後期課程1年　熊谷政仁​
東北大学工学部機械知能航空工学科　4年　深水一聖​
東北大学工学部機械知能航空工学科　3年　小野田誠

記事後半では発表では説明しきれなかった，詳しいQUBOの説明もしたいと思います．
かなり面白いテーマを，Domain Wall Encodingなどの最新技術で取り組んでいるので，読んでいただいているみなさんに何かいい気づきなどを与えることができれば嬉しいです．
では発表内容について紹介させてもらいます！



# 背景

## 環境保護のためCO2の削減
CO2削減目標が定められ，世界的に取り組まなければいけない課題となっています．

- パリ協定:気温上昇を2℃未満に,努力目標として1.5℃未満に抑える
- COP26:努力目標だった1.5℃が明確な目標に

また，先進国では2050年のカーボンニュートラル,2030年までの地球温暖化ガス排出量50%の削減が目標になっています．

## EVの世界的な普及・ガソリン車の廃止
2035年までに,世界的にガソリン車がEVに刷新される見通しであり，各国で政策レベルの議論が行われています．

- EU
    - 2035年までに,世界的にガソリン新車販売の禁止
    - ハイブリッド車も禁止
- 中国
    - 主要法規2021年にNEV規制が改訂
    - 2035年をめどに新車をすべて環境対応車
- アメリカ
    - カリフォルニア州で2035年までに州内の新車で内燃機関の廃
    - バイデン大統領が2030年までに国内販売車の新車40~50%電動化する大統領令を発令
- 日本
    - 2035年までに軽を含めて新車100%電動化方針を発表
## EVシェアの拡大
一方で，EVには次の2つの課題から，シェアを前提とした普及が叫ばれています．
- EVのバッテリーに使用するレアメタルの偏在や将来的な枯渇
    - 台数の確保が難しい
        - ガソリン車の4倍以上のレアメタルが必要
    - 価格が高くなりやすい
        - 300万を超える車体価格:補助金を含めてもガソリン車より高価
- シェアを前提としたEVの普及
    - 安価な値段で提供
    - 生産台数を抑え,消費資源や製造時のCO2削減
## 提供ソリューション
私たちは，東北大学主催の量子アニーリングソリューションコンテストにエントリーし，量子アニーリングを活用したEVシェアサービスとして，「EVerify ~EVeryone EVeryday EVerywhere~」に取り組みました．
![](https://i.imgur.com/zPWLaJy.jpg)


## 量子アニーリング×EVシェア
EVシェア実現のためには複雑な組み合わせ最適化問題を解く必要がありますが，量子アニーリングなら，最適化問題に特化した高速計算が期待できます．また，計算には冷却のためのごく僅かな電力しか使用しないので，環境にも優しい超省電力計算が期待できます．
量子アニーリングとEVシェアの組み合わせは，ハイパワー・クリーンエネルギーな最適化を実現に適していると考えています．
![](https://i.imgur.com/x93kUbg.jpg)



## 解決すべき課題
![](https://i.imgur.com/sld5Bk1.jpg)

### 住民目線の課題
![](https://i.imgur.com/0naUfhQ.jpg)

### 管理会社目線の課題
![](https://i.imgur.com/Guzynml.jpg)

### EV目線の課題
![](https://i.imgur.com/C89r8Ug.jpg)

# 最適化の全体像
最適化は，このワークフローの反復計算によって行われます．

- ユーザ：前日に，明日EVを利用したい時間を設定することで，競合が起きないようどのEVを使うべきか提案されます．
- 管理会社：ユーザ全体のスケジュール情報や，ユーザの居住地充電エリアを考慮して，明日どの充電エリアに何台のEVを配置しておけば良いか求めることができます．EVの必要台数を求めることは，設備投資にも役立ちます．
![](https://i.imgur.com/0gUpHKy.jpg)

**QUBO化については，次の流れで進めます．**
1. One Hot Encodingで定式化
2. One Hot Constraint部分を Domain Wall Constraintに変換
3. Domain Wall Constraint以外の部分のOne Hot Encodingで書かれた式をOne Hot Encodingに変換

One Hot Encoding ではなく，Domain Wall Encodingで最適式化をする理由は，One Hot Constraint部分のビット間相互作用を減らすためです．結合数の多いOne Hot Constraintはしばしば最適化の精度を低下させます．

## Domain Wall Encoding [[Chancellor, 2019](https://iopscience.iop.org/article/10.1088/2058-9565/ab33c2)]
- EVシェア最適化には多くの制約が必要
    - 道路状況の他に,充電時間や航続距離を考慮
        - 量子ビット間の相互作用が多くなる見込み
        - 量子アニーリングの精度が低下する可能性
    - Domain Wall Encodingによる離散変数表現
        - 従来広く用いられるOne Hot Encodingより少ない量子ビット数と相互作用数
![](https://i.imgur.com/fMwZ37r.jpg)
![](https://i.imgur.com/6CNsNdd.jpg)


| Performance       |One Hot Encoding  |Domain Wall Encoding|
|:----------------:|:-----------------:|:-----------:|
| 量子ビット数/スピン数             |    $$m$$           |   $$m-1$$     |
| 相互作用数         |            $$\frac{1}{2}m(m-1) $$ | $$m-2$$|
|　グラフ構造             |             完全グラフ |        線形グラフ |

Domain Wall Encodingにより,少ない量子ビット数と相互作用数で多制約最適化を実現できる可能性があります．




## EV割当・スケジューリング最適化
### 変数
まずは，One Hot Encodingによる定式化を行います．

- ユーザ$i$がEV$\alpha$を使用する$x_{i,\alpha}=1$
- ユーザ$i$がEV$\alpha$を使用しない$x_{i,\alpha}=0$

### 最適化に必要な要素
- 各ユーザのEV利用時間 $L_i$
- EV稼働可能時間 $T_{act}$
- 依存グラフ $E$ (あとで説明)

### 最適化の方針：EV利用時間Tαの均等化
いま，$T_{\alpha}=\sum_{\alpha}^{N−1} L_i x_{i, \alpha}$と EVの利用時間$T_{\alpha}$を定義します．

すると，EV利用時間の均等化は，次式によって行われます．
$
min \sum_0^{N−1}L_i x_{i, 0} \\
s.t.⁡T_{\alpha}　\leq　T_{0}
$

EV稼働可能時間 (充電時間)に関する制約は，次式で与えられます．
$
T_0 \leq T_{act}
$

### Log Encodingによる不等式制約のQUBO化
全てのEVの利用時間$T_{\alpha}$が$T_0$を超えない 

$$
T_{\alpha} \leq T_0
$$

$\overline{\Delta T(0,\alpha)}$を$T_{0}$と$T_{\alpha}$の差の上限とする



$$\therefore T_0−T_{\alpha}=\sum\limits_{i=0}^{N−1}L_i (x_{i,0}−x_{i,\alpha} ) \leq \overline{\Delta T(0,\alpha)} \because T_{\alpha}=\sum\limits_{0}^{N−1}L_i x_{i, \alpha} $$

Log Encoding により，

$0 \leq \sum_{i=0}^{{\lfloor\log_2 (\overline{\Delta T(0,\alpha)}-1 )\rfloor }}2^i𝑦_i \leq \overline{\Delta T(0,\alpha)}(𝑦_i\in{0, 1})$

を考える

$$
\overline{\Delta T(0,\alpha)} =  \sum\limits_{i=0}^{\lfloor\log_2 (\overline{\Delta T(0,\alpha)}-1 )\rfloor }2^i y_i +\sum\limits_{i=0}^{N−1}L_i(x_{i, 0}−x_{i, \alpha})
$$

等式制約なので，次の制約関数が得られる

$$
\min \Bigl\{⁡\sum_{\alpha=1}^{M−1}{\Delta T(0,\alpha)−\sum_{i=0}^{N−1}L_i (x_{i, 0}−x_{i, \alpha} )−\sum_{i=0}^{\lfloor\log_2 (\overline{\Delta T(0,\alpha)}-1 )\rfloor }2^i y_i}\Bigr\}^2
$$

$T_0$がEV稼働可能時間$T_{act}$を超えない$T_0　\leq T_{act}$についても，Log Encodingを用いた同様の式展開により，

$$
min\Bigl\{⁡(\Delta T(act,0)−T_{act}+\sum_{i=0}^{N−1}L_ix_{i, 0}−\sum_{i=0}^{\lfloor\log_2 (\overline{\Delta T(0,\alpha)}-1 )\rfloor }2^i z_i)^2 \Bigr\}
$$
ここまでで，充電時間を考慮したEV利用時間の均等化が完了しました．さらに，スケジュール・エリアの制約について考えます．


### エリア・スケジュールの同時依存回避
依存グラフ$E$について，次式を考えます．
$$
min⁡\sum\limits_{\alpha=0}^{M−1}\sum_{𝑖<𝑗}E_{i,j} x_{i, \alpha} x_{j, \alpha}
 $$

- エリア依存グラフ
エリアごとのユーザリスト：$U=\{0: (9, 4, 7, 8, 0), 1: (3, 5, 2), 2: (1, 6)\}$
行列$E_{area}$を考えて，ユーザ同士が別のエリアにいるなら1, 同じエリアなら0

- スケジュール依存グラフ
ユーザ同士がかぶっている利用時間数を要素にもつ行列$E_{schedule}$

２つの依存グラフを時間の単位に直して統合します．
$$E=T_{act} E_{area}+E_{schedule}$$

これにより，次にような最適化が行われます．
依存グラフの要素$E_{i,j}$が大きい時，$x_{i,\alpha}$と$x_{j, \alpha}$のどちらかは0 → スケジュールがかぶっている or エリアが異なる

$E_{i,j}=0$なら$x_{i,\alpha}$と$x_{j, \alpha}$が両方1になっても問題ない → スケジュールがかぶっていない and エリアが同じ


### 各ユーザをただ一つのEVに割り当てる制約
普通はOne Hot Constraintを使いますが，今回はDomain Wall Constraintを使います．これは相互作用数を削減するためです．
$$
 -\sum\limits_{i=0}^{N-1}\sum\limits_{\alpha=-1}^{M-2}(1-2b_{i,\alpha}-2b_{i,\alpha+1}+4b_{i,\alpha}b_{i,\alpha+1}) 
$$

### One Hot Encoding と Domain Wall Encoding の変換
$$x_{i, \alpha}=b_{i,\alpha−1}−b_{i, \alpha}$$
で行えます．


### QUBOの全体像
$$H_{objective}=\sum\limits_{0}^{N-1} L_i(b_{i,-1}-b_{i,0}) $$

$$H_{othercar}=\sum\limits_{\alpha=1}^{M-1} \Bigl\{\Delta T(0,\alpha)-\sum\limits_{i=0}^{N-1}L_i(b_{i,-1}-b_{i,0}-b_{i,\alpha-1}+b_{i,\alpha})-\sum\limits_{i=0}^{\lfloor\log_2 (\overline{\Delta T(0,\alpha)}-1 )\rfloor }2^iy_i \Bigr\}^2$$

$$H_{charging}= \biggl(\Delta T(act,0)-T_{act} + \sum\limits_{i=0}^{N-1}L_i(b_{i,-1}-b_{i,0}) -\sum\limits_{i=0}^{\lfloor\log_2 (\overline{\Delta T(act,0)}-1 )\rfloor }2^iz_i        \biggr) ^2$$

$$ H_{dependency} = \sum\limits_{\alpha=0}^{M-1}\sum_{i<j}E_{i,j}(b_{i,\alpha-1}-b_{i,\alpha})(b_{j,\alpha-1}-b_{j,\alpha})   $$

$$ H_{domain-wall}=-\sum\limits_{i=0}^{N-1}\sum\limits_{\alpha=-1}^{M-2}(1-2b_{i,\alpha}-2b_{i,\alpha+1}+4b_{i,\alpha}b_{i,\alpha+1}) $$
  
$$ H_{all}=AH_{objective}+B_1H_{other-car} + B_3H_{dependecy} + \kappa H_{damain-wall} $$ 

$$(A=1,B_1=1,B_2=1,\kappa = 24)$$

## 渋滞緩和最適化 [[Florian et.al., 2017](https://www.frontiersin.org/articles/10.3389/fict.2017.00029/full)]
Volks Wagenの有名な論文と同じQUBOを使います．ある車$i$に対して$k$個の経路パターンを与え$x_{i,k}$を選ぶときは$x_{i,k} = 1$、選ばないときは$x_{i,k}=0$として考えます．

コスト関数は、ある経路$(i,k)$について、道路$e$を通るかどうかを示す$C_{e, (i,k)}$を用意します．この$C_{e, (i,k)}$は$0$または$1$を持ち、道路$e$を通る場合は$C_{e, (i,k)}=1$、道路$e$を通らない場合は$C_{e, (i,k)}=0$とします．


One Hot Encoding を用いたQUBOは次にようになります．
 $$ H= \sum_e  \biggl( \sum\limits_{i=0}^{N-1}\sum\limits_{k=0}^{K-1}C_{e(i,k)}x_{i,k} \biggr)^2+\lambda\sum\limits_{i=0}^{N-1}\biggl( \sum\limits_{k=0}^{K-1}x_{i,k-1} \biggr)^2$$


これにDomain-Wall　Encodingを適用すると,One-Hot ConstraintがDomain-Wall Constraintになります．さらに，$x_{i,k}=b_{i,k-1}-b_{i,k}$を適用します．
$E(x) = \sum_{e} \left( \sum_{i,k} C_{e, (i,k)}(b_{i,k-1}-b_{i,k}) \right)^2 - \kappa \sum_{i=0}^{N-1}  \sum_{k=-1}^{K-2} (1-2b_{i,k})(1-2b_{i,k-1})
$

  

# 仙台市泉パークタウンでのデモンストレーション
- 泉パークタウンのスマートシティ計画 ([参考資料](https://www.mlit.go.jp/toshi/tosiko/content/001397752.pdf))
    - 住民目線による都市OS最適化を目指し,持続可能な街全体のアーキテクチャを構築
    - 交通・物流面では,公共交通機関の効率化とEVシェア拠点の構築
- デモンストレーション
    - 泉パークタウンの住民を想定してユーザを選択し，EV割り当て・スケジューリング最適化
    - 実際の地図上で，渋滞緩和最適化を行い，交通量を可視化

コードは[Colab](https://colab.research.google.com/drive/1KUv8FcYkjLqi-6LpwpWFZVng5zyczQbR?usp=sharing
)で公開中です．ぜひ試してみてください！


## EV割り当て・スケジューリング最適化
![](https://i.imgur.com/d7SUdDS.jpg)
![](https://i.imgur.com/SOLWijd.jpg)

## 渋滞緩和最適化
![](https://i.imgur.com/O19jCOX.jpg)
![](https://i.imgur.com/vw5XKIt.jpg)





# まとめ
この記事ではソリューションコンテストで発表した内容を，発表動画には入りきらなかった部分の説明も込めて詳細に紹介してきました．
この記事を書いている2021/12/10現在ではまだ結果は出ていませんが，結構いい線行ってるのではないかと思います．
結果がわかれば情報をアップデートできればいいなと思っております．

また，繰り返しですがこの記事の内容について，[YouTube動画](https://www.youtube.com/watch?v=6er7f41l-0g&ab_channel=TohokuForumforCreativity)が出ています！是非ご覧ください！

最後にQuantumPredatorのメンバーと，ここまでこの記事を読んでいただいた皆さん，本当にありがとうございました．



量子アニーリングを用いたEVシェアの最適化：量子アニーリングソリューションコンテストより

EVの世界的な普及・ガソリン車の廃止

EV割当・スケジューリング最適化

最適化の方針：EV利用時間Tαの均等化

Log Encodingによる不等式制約のQUBO化

エリア・スケジュールの同時依存回避

各ユーザをただ一つのEVに割り当てる制約

One Hot Encoding と Domain Wall Encoding の変換

渋滞緩和最適化 [Florian et.al., 2017]

仙台市泉パークタウンでのデモンストレーション

EV割り当て・スケジューリング最適化

グラフ彩色問題 (PyQUBOを用いたQUBO生成)

アニーリングマシンで機械学習

## 見出し
```
# 見出し1
## 見出し2
### 見出し3
```

## 文字装飾
```
**太字**
*斜体*
~~打ち消し線~~
`コード`
$ax^2+bx+c$
```
**太字**
*斜体*
~~打ち消し線~~
`コード`
$ax^2+bx+c$

## リンク
```
[リンクテキスト](URL)
```
[ANCAR](https://ancar.app)

## 画像
```
![代替テキスト](画像URL)
```
![](https://firebasestorage.googleapis.com/v0/b/ancar-8c06d.appspot.com/o/articles%2FU10-4266-6317-93%2F17bb37fd1ef1eb?alt=media&token=a2a08166-c2d4-4689-9e92-2bf55d277d97=250x)


## リスト
```
// 順序なし
- 項目A
- 項目B
- 項目C
```
- 項目A
- 項目B
- 項目C


```
// 順序なし
1. 項目1
2. 項目2
3. 項目3
```
1. 項目1
2. 項目2
3. 項目3


## テーブル
```
| Head | Head | Head |
| ---- | ---- | ---- |
| Text | Text | Text |
| Text | Text | Text |
```

| Head | Head | Head |
| ---- | ---- | ---- |
| Text | Text | Text |
| Text | Text | Text |


## コード
> \```
> print("Hello World")
> \```

```
print("Hello World")
```

## 数式
```
$$
x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}
$$
```

$$
x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}
$$

## 引用
```
> 引用文
```
> 引用文

## 水平線
```
---
```

---

投稿コンテンツ（新着）