組合せ最適化でチーム分けする（焼きなまし法）

## はじめに
[前回の記事](http://qiita.com/matsulib/items/898873b73d584c7dcb8b)でチーム分けの問題を平均偏差最小化問題として定式化し、Python+PuLPで厳密に解く方法を紹介しました。その中で、組合せ最適化問題は問題が大規模になると厳密に解くことが難しくなることを述べました。

そこで、今回はPython＋simannealを使った焼きなまし法による近似解法を紹介します。

## simannealとは

simannealは焼きなまし法のためのパッケージです。

**インストール方法**

```bash
pip install simanneal
```

**基本的な使い方**

* Annealerクラスのサブクラスを作成する。
    * move()をオーバーライドし、探索点の移動ルールを実装
    * energy()をオーバーライドし、目的関数を実装
* 作成したクラスのコンストラクタに初期解を与えて問題オブジェクトを作成する。
* 問題オブジェクト.anneal()で焼きなまし法を実行する。
* 問題オブジェクト.stateで解を取得する。

最適化問題の制約はmove()とenergy()の中で表現する必要があります。

後述の「例題：Python＋simannealによる最適化」も参照して下さい。

## 解きたい問題
前回と同じ問題を扱います。
背景は前回の記事を参照してもらうとして、以下のような組合せ最適化問題を解きます。

```math
\begin{align}
\text{minimize} \quad &\sum_{j=1}^m \Bigl(\sum_{k=1}^l \Bigl|\sum_{i=1}^n a_{i,k} x_{i,j} - \mu_j \Bigr|\Bigr) + 1.5\sum_{j=1}^m\Bigl|\sum_{i=1}^n b_i x_{i,j} -\nu\Bigr| \tag{2.1}\\
\text{subject to} \quad &\mu_j = \frac{1}{l}\sum_{k=1}^l\sum_{i=1}^n a_{i,k}x_{i,j} &j=1,\ldots,m \tag{2.2}\\
& \nu = \frac{1}{m} \sum_{j=1}^m \sum_{i=1}^n b_i x_{i,j} \tag{2.3}\\
&\sum_{j=1}^m x_{i,j} = 1 &i=1,\ldots,n \tag{2.4}\\
& x_{i,j} \in \{0,1\} &i=1,\ldots,n;\, j=1,\ldots,m \tag{2.5}\\
\end{align}
```
前回はここから絶対値を外して線形計画問題にするような変形を行いましたが、焼きなまし法ではそのままで大丈夫です。

## 例題：Python＋simannealによる最適化

前回と同じ例題を扱います。

* 6人を3チームに分けます。
* コミュニケーション能力の種類は4種類です。

```python:grouping.py
from simanneal import Annealer
import random

class GroupingProblem(Annealer):
    
    def __init__(self, init_state):
        super(GroupingProblem, self).__init__(init_state)  # important!

    # 探索点の移動ルール
    def move(self):
        # ランダムに選択したメンバーm1をチームt1からチームt2に移す
        m1 = random.choice(list(range(nmembers)))
        t1 = self.state[m1].index(1)
        t2 = random.choice(list(range(nteams)))
        self.state[m1][t1], self.state[m1][t2] = self.state[m1][t2], self.state[m1][t1]

    # 目的関数
    def energy(self):
        v = 0
        nu = sum(sum(b[i] * self.state[i][j] for i in range(nmembers)) for j in range(nteams)) / nteams
        for j in range(nteams):
            mu_j = sum(sum(scores[i][k] * self.state[i][j] for i in range(nmembers)) for k in range(nskills)) / nskills
            for k in range(nskills):
                v += abs(sum(scores[i][k] * self.state[i][j] for i in range(nmembers)) - mu_j)
            v += 1.5 * abs(sum(b[i] * self.state[i][j] for i in range(nmembers)) - nu)
        return v


if __name__ == '__main__':
    # 例題
    teams = ['A', 'B', 'C']
    members = ['P', 'Q', 'R', 'S', 'T', 'U']
    skills = ['Controller', 'Promoter', 'Supporter', 'Analyzer']
    scores = [[6, 0, 1, 3],
              [2, -1, 3, 5],
              [2, 4, 0, 0],
              [3, 4, 5, 0],
              [0, 2, 1, 4],
              [2, 3, -1, 1]]

    nteams = len(teams)  # チーム数
    nmembers = len(members)  # メンバー数
    nskills = len(skills)  # 能力種別数
    b = [sum(ai) for ai in scores]  # 各メンバーのスキルの総和

    # 初期割当て
    init_state = [[0 for j in range(nteams)] for i in range(nmembers)]
    for i in range(nmembers):
        init_state[i][0] = 1  # 最初は全員Aチームに所属
    
    prob = GroupingProblem(init_state)
    prob.steps = 100000
    prob.copy_strategy = "deepcopy"
    prob.anneal()  # 焼きなまし実行

    for i,s in enumerate(prob.state):
        print(members[i], teams[s.index(1)])
```
これを実行すると焼きなましの経過が表示されます。
![3.gif](https://qiita-image-store.s3.amazonaws.com/0/141719/4c2263f3-c9de-1660-8596-386d8e43cea4.gif)

前回と同じ解が得られました。

```bash
$ pypy grouping.py
 Temperature        Energy    Accept   Improve     Elapsed   Remaining
     2.50000         23.50    35.60%     0.70%     0:00:04     0:00:00
('P', 'A')
('Q', 'C')
('R', 'C')
('S', 'B')
('T', 'A')
('U', 'B')
```
simanneal が表示しているそれぞれの項目について補足します。

* Temperature: 温度。この値が大きいときほど、状態変化が激しくなる。
* Energy: エネルギー。目的関数の値。
* Accept: 遷移確率。この値が大きいと悪い状態へも遷移する（局所解から脱出するため）。この値はエネルギーと温度に依存する。
* Improve: 遷移した際の目的関数の改善率。
* Elapsed: 経過時間。
* Remaining: 残り時間。

最終的にはTemperatureとImproveが小さくなり解が収束します。

## おわりに
組合せ最適化問題をPython＋simannealで解く方法を紹介しました。
simannealは簡単に焼きなまし法を実行できるとても便利なパッケージです。
また、simannealはpypyでも実行できます。イテレーションごとの操作(move()やenargy()の中身)が単純なものならばpypyを利用することで高速化が見込めます。

注：pypyでnumpy.arrayを使うと遅くなります。pypyを使う場合はnumpy.arrayのベクトル演算を使うよりリストでmapやfor文を回したほうが速いです。

###参考

* [perrygeo/simanneal: Python module for Simulated Annealing optimization](https://github.com/perrygeo/simanneal)
* [焼きなまし法 - Wikipedia](https://ja.wikipedia.org/wiki/%E7%84%BC%E3%81%8D%E3%81%AA%E3%81%BE%E3%81%97%E6%B3%95)
* [Numpy を使わない多次元配列計算](http://qiita.com/ponnhide/items/c919f3bc549d1228c800)