数理最適化で、確率に応じて委員を決めよう

## 問題

クラスの委員が5人だけ決まっていません（体育が2人、図書が1人、保健が2人）。
時間もないので抽選することになりました。
なるべく生徒の**希望に応じた確率**で決めてください。

### 希望

* 生徒1：体育か図書を希望。第1希望が体育
* 生徒2：体育か保健を希望。第1希望が体育
* 生徒3：保健か体育を希望。第1希望が保健
* 生徒4：図書か保健を希望。第1希望が図書
* 生徒5：保健か体育を希望。第1希望が保健

ここでは、**第1希望の確率を2/3、第2希望の方を1/3**としましょう。
また、前提として、委員ごとの確率の和は委員の人数に一致することとします。

## 考え方

この割当問題を数理最適化でモデル化します。
まず、確率から何らかの重みを算出します。
目的関数は、その重みを上限とする一様乱数を係数とします。
目的関数を最大化することで、確率で割り当たるようにしてみましょう。

### 重みと確率

ここから少し寄り道をします。

2つの重みを2と3としましょう。
この重みから2つの一様乱数（範囲は0〜2と0〜3）aとbを生成します。
aとbを2次元で表すと、取りうる範囲は2 x 3で面積は6です。
aが大きくなる面積は2で、bが大きくなる面積は4です。
この2と4は、次のように関数`func2()`で計算できます。

**In**

```python
import numpy as np

def func2(a, b):
    p = a * a / 2
    yield p
    q = p + a * (b - a)
    yield q

weight = np.array([2, 3])
print(*func2(*weight))
```

**Out**

```
2.0 4.0
```

実際に60万組の乱数を発生して確認してみましょう。

**In**

```python
n = 600_000
rng = np.random.default_rng(0)
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
```

**Out**

```
(array([0, 1]), array([199555, 400445]))
```

おおよそ2対4になっています。

関数`func2()`の逆関数`r_func2()`は、式変形することで次のように作成できます。

**In**

```python
def r_func2(p, q):
    a = (p * 2)**(1 / 2)
    yield a
    b = (p + q) / a
    yield b

rate = np.array([2, 4])
print(*(r_func2(*rate)))
```

**Out**

```
2.0 3.0
```

`weight`と同じ値が逆算されました。

重みの個数が3つのときは、次のようになります。これは、3次元の体積から計算できます。

**In**

```python
def func3(a, b, c):
    p = 0 + 1 * (a**3 - 0**3) / 3
    yield p
    q = p + a * (b**2 - a**2) / 2
    yield q
    r = q + a * b * (c**1 - b**1) / 1
    yield r

def r_func3(p, q, r):
    a = (3 * (p - 0) + 0**3)**(1 / 3)
    yield a
    b = (2 * (q - p) / a + a**2)**(1 / 2)
    yield b
    c = (1 * (r - q) / a / b + b**1)**(1 / 1)
    yield c

rate = np.array([1, 2, 3])
print(*r_func3(*rate))
```

**Out**

```
1.4422495703074083 1.8619361889584654 2.2343234267501586
```

逆関数が正しいことを確認しましょう。

**In**

```python
weight = np.array([*r_func3(*rate)])
print(np.array([*func3(*weight)]))
```

**Out**

```
[1. 2. 3.]
```

実際に60万組の乱数を発生して確認してみましょう。指定した比率になっています。

**In**

```python
n = 600_000
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
```

**Out**

```
(array([0, 1, 2]), array([100019, 199960, 300021]))
```

重みの個数が4つのときは、次のようになります。計算式は、なんとなく規則性がわかるかと思います。

**In**

```python
def func4(a, b, c, d):
    p = 0 + 1 * (a**4 - 0**4) / 4
    yield p
    q = p + a * (b**3 - a**3) / 3
    yield q
    r = q + a * b * (c**2 - b**2) / 2
    yield r
    s = r + a * b * c * (d**1 - c**1) / 1
    yield s

def r_func4(p, q, r, s):
    a = (4 * (p - 0) / 1 + 0**4)**(1 / 4)
    yield a
    b = (3 * (q - p) / a + a**3)**(1 / 3)
    yield b
    c = (2 * (r - q) / a / b + b**2)**(1 / 2)
    yield c
    d = (1 * (s - r) / a / b / c + c**1)**(1 / 1)
    yield d

rate = np.array([1, 2, 3, 4])
weight = np.array([*r_func4(*rate)])
n = 1_000_000
rnd = rng.random((n, weight.size))
print(np.unique((rnd * weight).argmax(axis=1), return_counts=True))
```

**Out**

```
(array([0, 1, 2, 3]), array([100125, 200424, 299249, 400202]))
```

## 解答

各生徒の希望確率を次のようにします。

```python
wish_probs = [
    {"体育": 2 / 3, "図書": 1 / 3},
    {"体育": 2 / 3, "保健": 1 / 3},
    {"保健": 2 / 3, "体育": 1 / 3},
    {"図書": 2 / 3, "保健": 1 / 3},
    {"保健": 2 / 3, "体育": 1 / 3},
]
```

この確率の割合に対する重みは、次のように`2.0`と`1.5`です。

**In**

```python
print(*r_func2(2, 1))
```

**Out**

```
2.0 1.5
```

ここからは、実際に問題を解いていきます。

最初に、モデル作成で使う関数などを定義します。

```python
from dataclasses import make_dataclass
import polars as pl
import numpy as np
from mip import Model, maximize, xsum


def xdot(a, b):
    """内積"""
    a = a.to_numpy()
    b = b.to_numpy()
    return xsum(a * b)


def int_value(v):
    """整数解"""
    return v.to_numpy().astype(float).round().astype(int)


rng = np.random.default_rng(0)
```

次に、数理モデルを作成します。

```python
m = Model()
df = pl.from_dicts(
    {"生徒": i, "委員": committee, "確率": prob}
    for i, wish_prob in enumerate(wish_probs)
    for committee, prob in wish_prob.items()
)
df = df.with_columns(
    係数=pl.when(pl.col("確率") > 0.5).then(2.0).otherwise(1.5),
    Var=m.add_var_tensor((len(df),), "Var", var_type="B"),
    Val=0,
)
col = make_dataclass("Col", df.columns)(*map(pl.col, df.columns))

for _, group in df.group_by(col.委員):
    # 委員ごとの人数
    m += xsum(group["Var"]) == group["確率"].sum()
for _, group in df.group_by(col.生徒):
    # 生徒ごとに1つ選択
    m += xsum(group["Var"]) == 1

m.verbose = 0
```

最後に、目的関数を変えながら1000回解いて、生徒ごとの割合を見てみましょう。

**In**

```python
result = []
for _ in range(1000):
    coef = df["係数"] * rng.random(len(df))
    m.objective = maximize(xdot(coef, df["Var"]))
    m.optimize()
    df = df.with_columns(Val=int_value(df["Var"]))
    result.append(df.filter(col.Val == 1).sort("生徒")["委員"].to_list())
for lst in zip(*result):
    print(np.unique(lst, return_counts=True))
```

**Out**

```
(array(['体育', '図書'], dtype='<U2'), array([699, 301]))
(array(['体育', '保健'], dtype='<U2'), array([589, 411]))
(array(['体育', '保健'], dtype='<U2'), array([336, 664]))
(array(['保健', '図書'], dtype='<U2'), array([301, 699]))
(array(['体育', '保健'], dtype='<U2'), array([376, 624]))
```

生徒同士は独立ではないので、必ずしも指定した割合になるとは限りませんが、近い値にはなっているようです。

## まとめ

確率で割り当てる問題を数理最適化で解いてみました。
各委員の人数は決まっているため、生徒ごとに独立に決めることはできません。
数理モデルを使うことで、なるべく確率ごとに選んだり、委員の人数を指定できました。

以上