組合せ最適化でチーム分けする（平均偏差最小化）

# はじめに
[組合せ最適化でチーム分けする](http://qiita.com/SaitoTsutomu/items/f4478dfbc3c1cf6425e3) で紹介されている問題を別のモデルで解いてみます。

この問題は、バラツキを最小化する組合せ最適化問題です。
バラツキ具合の指標としては一般的に分散がよく使われますが、分散は二次関数で表されるため一次関数を最小化する線形計画問題に比べると難しくなってしまいます。そこで元の記事では最大値と最小値の差（範囲）を最小化する線形計画問題として定式化しています。

一方、範囲の代わりに平均偏差を最小化する方法もあります。
そこで本記事では、まず範囲と平均偏差のモデルの違いを比較します。
最後に平均偏差最小化問題の例題をPython＋PuLPで解いてみます。

# 問題設定

* 前提条件
 * $n$人のメンバーを$m$個のチームに分ける。
 * 各メンバーが持つコミュニケーション能力の種類は$l$個あり、それぞれ点数が振られている。
* 目標
 * チーム内で能力が偏らないようにする。
 * チーム間の能力差も小さくする。

したがって、目的関数の基本方針は以下のようになります。

$$
\text{minimize} \quad \sum_{j=1}^m (\text{チーム$j$内のバラツキ}) + (\text{チーム間のバラツキ}) \qquad \tag{0.1}
$$

また、どちらかのバラツキに重みを付けてバランスを取ることもできます。

# 定式化
問題を定式化していきます。
まず、以下の変数を定義します。

```math
\begin{align}
n &: メンバー数 \\
m &: チーム数 \\
l &: スキル数 \\
x_{i,j} &: メンバーiのチームjへの割当 (\in \{0,1\}) \\
a_{i,k} &: メンバーiのスキルkの点数 \\
b_i &: メンバーiの全スキルの総和 (=\sum_{k=1}^l a_{i,k}) \\
\end{align}
```

## 範囲最小化問題

バラツキの指標として、最大値と最小値の差（範囲）を用いて、それを最小化します。

ここで、

```math
\begin{align}
y_{j,\min}, y_{j,\max} &: チームj内のスキルの最小と最大 \\
z_\min, z_\max &: チーム間の最小と最大
\end{align}
```

とすると、以下のように範囲最小化問題として定式化できます。

```math
\begin{align}
\text{minimize} \quad &\sum_{j=1}^m (y_{j \max} -y_{j \min}) + 1.5 (z_\max - z_\min) \tag{1.1}\\
\text{subject to} \quad & y_{j \min} \leq \sum_{i=1}^n a_{i,k} x_{i,j} \leq y_{j \max} &j=1,\ldots,m;\, k=1,\ldots,l \tag{1.2}\\
& z_\min \leq \sum_{i=1}^n b_i x_{i,j} \leq z_\max &j=1,\ldots,m \tag{1.3}\\
&\sum_{j=1}^m x_{i,j} = 1 &i=1,\ldots,n \tag{1.4}\\
& x_{i,j} \in \{0,1\} &i=1,\ldots,n;\, j=1,\ldots,m \tag{1.5}\\
\end{align}
```
式(1.4)はどのメンバーもどれか一つのチームに所属することを表す制約式です。
元の記事にあわせて、チーム間のバラツキの方に重み1.5を掛けています。
変数の数も少なく、とてもシンプルに定式化できます。

## 平均偏差最小化問題

$x_i$を個々の値、$ \bar{x} $をその平均とすると、分散と平均偏差は次のように表せます。

* $ 分散=\frac{1}{n}\sum_i^n (x_n - \bar{x})^2 $ 
* $ 平均偏差=\frac{1}{n}\sum_i^n |x_n - \bar{x}| $

平均偏差には絶対値が含まれており微分不可能で計算が面倒などの理由で避けられていますが、最適化問題の目的関数に絶対値が含まれる場合は簡単に外すことができます。

平均偏差を使って式(0.1)を書き直すと、以下のようになります。

```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}
```

このような絶対値最小化問題は補助変数$y_{j,k}, \,z_j$を導入することで、以下のように線形計画問題に変形することができます。

```math
\begin{align}
\text{minimize} \quad &\sum_{j=1}^m \Bigl(\sum_{k=1}^l y_{j,k} + 1.5 z_j\Bigr) \tag{2.6}\\
\text{subject to} \quad & -y_{j,k} \leq \sum_{i=1}^n a_{i,k} x_{i,j} - \mu_j \leq y_{j,k} & j=1,\dots,m;\,k=1,\ldots,l \tag{2.7}\\
& -z_j \leq \sum_{i=1}^n b_i x_{i,j} -\nu \leq z_j & j=1,\ldots,m \tag{2.8}\\
&\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.9}\\
& \nu = \frac{1}{m} \sum_{j=1}^m \sum_{i=1}^n b_i x_{i,j} \tag{2.10}\\
&\sum_{j=1}^m x_{i,j} = 1 &i=1,\ldots,n \tag{2.11}\\
& x_{i,j} \in \{0,1\} &i=1,\ldots,n;\, j=1,\ldots,m \tag{2.12}\\
\end{align}
```
範囲最小化のモデルと比較すると、変数の数が多くなり、若干煩雑なモデルになりました。

# 例題：Python＋PuLPによる最適化
PythonのPuLPというパッケージを使って例題を解いてみます。

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

PuLPではソルバーを選択できますが、ここではデフォルトのCoin-CBCを使います。

```python
import numpy as np
import pandas as pd
from pulp import *
from ortoolpy import addvar, addvars, addbinvars

# 例題
teams = ['A', 'B', 'C']
members = ['P', 'Q', 'R', 'S', 'T', 'U']
skills = ['Controller', 'Promoter', 'Supporter', 'Analyzer']
scores = pd.DataFrame([[6, 0, 1, 3],
                       [2, -1, 3, 5],
                       [2, 4, 0, 0],
                       [3, 4, 5, 0],
                       [0, 2, 1, 4],
                       [2, 3, -1, 1]],
                      index=members,
                      columns=skills)

nteams = len(teams)  # チーム数
nmembers = len(scores.index)  # メンバー数
nskills = len(scores.columns)  # 能力種別数

print(scores)

#    Controller  Promoter  Supporter  Analyzer
# P           6         0          1         3
# Q           2        -1          3         5
# R           2         4          0         0
# S           3         4          5         0
# T           0         2          1         4
# U           2         3         -1         1
```

## 範囲最小化
[元の記事](http://qiita.com/SaitoTsutomu/items/f4478dfbc3c1cf6425e3)を参照して下さい。

## 平均偏差最小化

```python
# 平均偏差最小化

m = LpProblem()  # 数理モデル
x = pd.DataFrame(addbinvars(nmembers, nteams), index=members, columns=teams)  # 割当
y = pd.DataFrame(addvars(nteams, nskills), index=teams, columns=skills)  # チーム内の平均偏差
mu = pd.DataFrame(addvars(nteams), index=teams)   # チーム内の平均
z = pd.DataFrame(addvars(nteams), index=teams)  # チームごとの平均偏差
nu = addvar()  # 全チームの平均

m.setObjective(lpSum([lpSum(y.loc[j]) + 1.5*z.loc[j] for j in teams]))  # 目的関数

m.addConstraint(lpSum(np.dot(scores.sum(1), x)) / nteams == nu)
for j in teams:
    m.addConstraint(lpDot(scores.sum(1), x[j]) - nu <= z.loc[j])
    m.addConstraint(lpDot(scores.sum(1), x[j]) - nu >= -z.loc[j])
    m.addConstraint(lpSum(np.dot(x[j], scores)) / nskills == mu.loc[j])
    for k in skills:
        m.addConstraint(lpDot(scores[k], x[j]) - mu.loc[j] <= y.loc[j,k])
        m.addConstraint(lpDot(scores[k], x[j]) - mu.loc[j] >= -y.loc[j,k])
for i in members:
    m.addConstraint(lpSum(x.loc[i]) == 1)  # どこかのチームに所属
        
m.solve()  # 求解

if m.status == 1:
    for i, xi in enumerate(np.vectorize(value)(x).tolist()):
        print((members[i], teams[xi.index(1)]))
else:
    print('最適解を得られなかった。')

# 結果：(メンバー, チーム)
# ('P', 'B')
# ('Q', 'A')
# ('R', 'A')
# ('S', 'C')
# ('T', 'B')
# ('U', 'C')
```
元の記事と同じ解が得られました。
別のモデルを解いているので同じになったのはたまたまです。

# おわりに
平均偏差最小化問題を取り上げました。
このモデルは、例えばポートフォリオ最適化で使えます。ポートフォリオ最適化では期待リターンのバラツキをリスクとみなし、リスク最小化問題として定式化しますが、その際に平均偏差を最小化することがあります。
ただし決定変数に整数が含まれる整数計画問題（組合せ最適化問題、離散最適化問題）の場合は、問題の規模が大きくなるとあっという間に解けなくなります。
今回取り上げた問題も「20人を5チームに分ける」くらいになっただけで僕のPCでは現実的な時間で解くのが厳しくなりました（一晩放置しても解けない）。ソルバーとしてGurobiやCPLEXのような商用ソルバーを利用することで特にマルチコア環境でかなりの改善が見込めますが、それでも大規模な問題を解くのは難しいでしょう。その場合はモデルを見直すか、近似解法を使うことが考えられます。

【続き】今回の問題の近似解法として、焼きなまし法のためのPythonパッケージsimannealを使って解く方法を紹介します。
[組合せ最適化でチーム分けする（焼きなまし法）](http://qiita.com/matsulib/items/bd50af2e2bc1e48522cd)

### 参考

* [Python+PuLPによるタダで仕事に使える数理最適化](http://qiita.com/samuelladoco/items/703bf78ea66e8369c455)
* [組合せ最適化を使おう](http://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721)