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

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

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で自動化する

# はじめに
初めまして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」に関する投稿コンテンツ

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

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

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