最適かつ公平に資源を配分したい - ナッシュ交渉解の解釈とLPソルバーを用いた計算

　この記事は [数理最適化 Advent Calendar 2024](https://qiita.com/advent-calendar/2024/mathematical-optimization) の 12 日目の記事です．
　12月12日は国際中立デー(International Day of Neutrality)だそうですね．特に関係はありませんが，資源配分の公平性について紹介したいと思います．

# 資源配分の問題設定

　需要家に対して資源を配分したいという状況を想定します．資源には複数の種類が存在し，それぞれの供給量が決められているとします．需要家も複数存在し，それぞれが配分された資源の量に対しての利得関数を持っており，利得関数は既知である（需要家が正直に申告した）ものとします．
　このときに，各需要家の利得を最大化し，かつ公平な配分量を決定する問題を考えます．

　このような問題の例としては，以下のようなものが考えられます．

* 店舗業務の従業員への割り当て
    * 店舗では，陳列・レジ・清掃など複数の業務があり，従業員ごとに担当したい業務が異なる
        * 力仕事は苦手・接客がやりたい・きれい好きなど
    * なるべく従業員の要望を叶えたいが，みんながやりたがっている業務が特定の従業員に集中するといった不公平感がないようにしたい
* 広告機会(TVCM枠・リスティング広告など)の広告主への配分
    * 広告機会ごとに視聴するユーザーが異なり，広告主ごとにターゲットとするユーザーが異なる
    * 広告主に対して，ターゲットとしているユーザーの視聴率が高い広告を配分することで効率よく利得を大きくすることができるが，特定の広告主の利得が小さくなってしまうといった不公平がないようにしたい

　本記事では，資源の需要家への配分量は連続量であり，需要家の利得関数が配分量に対して線形である場合に限定して，以下の問題を扱うことにします．

```math
\begin{align}
\text{max} \ & \ y_1, \cdots , y_n \\
\text{s.t.} \ & \ \sum_{i=1}^m w_{ij} x_{ij} = y_j \\
& \ \sum_{j=1}^n x_{ij} \leq s_i \\
& \ x_{ij} \geq 0 \\
& \ y_j \in \mathbb{R}
\end{align}
```

　資源を$i \in \lbrace 1, \cdots, m \rbrace$，需要家を$j \in \lbrace 1, \cdots, n \rbrace$とし，資源$i$の需要家$j$への配分量を変数$x_{ij}$，需要家$j$の利得を変数$y_j$で表すことにします．$s_i \geq $は資源$i$の供給量，$w_{ij} \geq 0$は配分量$x_{ij}$の利得への寄与(需要家$j$の利得関数は$f_j(x)=\sum_i^m w_{ij} x_{ij}$)として与えられているとします．
 
　この問題は多目的最適化になっており，このままでは配分量$x$を決定することができないため，以下ではこの問題を解ける構造に変換するための 2 つの方法を考えて，次の具体的な例を使って結果を確認していきます．

**【具体例】**
　3 つの資源カレーライス・ラーメン・プロテインがそれぞれ 1 食分存在し，それらの資源を 2 人の需要家エネルギッシュ・マッスルに配分する状況を考えます．エネルギッシュの利得は配分された資源のエネルギーの合計で，マッスルの利得は配分された資源のたんぱく質の合計であるとします．また，需要家ごとの支払額が与えられているとします．（資源が単品メニューなので忘れそうになりますが，ここでも資源の配分量は整数でなくても良いとしています）

![image1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/8ca879fc-e44d-a03f-bde7-c01ebe0e96d1.png)

　2 人の需要家は全く異なる価値観で資源を評価しているため，適切な配分方法は自明ではありません．

# 資源配分問題を解く

## 利得の線形和を目的関数とする方法

　複数の指標を同時に扱いたいときには，各指標に対して適当な重み $c$ を掛けてそれらの和を目的関数とする方法がよく使われます．

```math
\begin{align}
\text{max} \ & \ \sum_{j=1}^{n} c_j y_j \\
\text{s.t.} \ & \ \sum_{i=1}^m w_{ij} x_{ij} = y_j \\
& \ \sum_{j=1}^n x_{ij} \leq s_i \\
& \ x_{ij} \geq 0 \\
& \ y_j \in \mathbb{R}
\end{align}
```

　上に挙げた具体例では，利得の単位を揃えつつ支払額も考慮した重みとして，例えば以下が考えられます．

```math
\begin{align}
c_\text{エネルギッシュ}&=\frac{1000}{600 + 500 + 100}\\
c_\text{マッスル}&=\frac{1000}{10 + 20 + 20}
\end{align}
```

　この重みで問題を解くと，カレーライスとラーメンは全てエネルギッシュに配分され利得は 1100 kcal，プロテインは全てマッスルに配分されて利得は 20g という結果が得られます．

![graph1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/3bf9bf11-8594-0a47-1c40-bac8f8a855f4.png)

　パレート最適解なので利得の最大化は達成できていますが，公平な配分になっていると言えるでしょうか？ 
　公平とは何かという統一された定義は存在しませんが，本記事では，需要家から見た資源の単価(需要家にとってのコストパフォーマンスの逆数)を使って公平性を考えてみたいと思います．

　この結果では，エネルギッシュは ￥1000 を支払って 1100 kcal を得たので，エネルギッシュにとってのカレーライスの単価(1食分の価格)は
 ```math
 600 \text{kcal} / 1100 \text{kcal} \times \text{￥}1000 = \text{￥}545 \ \ \text{(小数点以下は四捨五入)}
 ```
であるとみなすことができます．同様にして，各需要家にとっての各資源の単価を計算すると以下のようになります．

![image2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/a4fdbf8e-c475-5f5d-ead3-e6756bab45d1.png)

　ここでラーメンの単価に注目すると，エネルギッシュにとっては ￥455，マッスルにとっては ￥1000 となっており，「マッスルの方がエネルギッシュよりもラーメンを高く評価しているにも関わらず，ラーメンは全てエネルギッシュに配分されている」という状況が発生しています．

　この状況は，需要家間での転売が成立し得るという意味で公平ではなく，もしもエネルギッシュが配分されたラーメン 1 食分を ￥999 でマッスルに転売すると，エネルギッシュは実質タダ同然(￥1)でカレーライス 1 食分を得られることになってしまいます．
 
 ![image3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/cadf381a-2b41-98e9-e56e-38fb84accbc3.png)

<details><summary><b>【参考】マッスルの支払額を￥1100としたケース</b></summary>

　マッスルの支払額を ￥1100 に増やして目的関数の重みを以下のように変更してみましょう．

 ```math
\begin{align}
c_\text{エネルギッシュ}&=\frac{1000}{600 + 500 + 100}\\
c_\text{マッスル}&=\frac{\color{red}{1100}}{10 + 20 + 20}
\end{align}
```

　この重みでは，カレーライス1食分がエネルギッシュに配分され，ラーメン1食分とプロテイン1食分とがマッスルに配分される結果が得られます．

![graph2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/a479b6ea-bbf9-9f45-aadd-bdf423d35c08.png)

　マッスルの支払額を10%増やしただけですが，それぞれの利得は大きく変動し，先ほどとは逆にマッスルからエネルギッシュへの転売が成立しています．

![image4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/56a2dfa0-6f0f-54d7-630c-3998ca371ca0.png)

</details>

　以上のように，需要家の利得の線形和を目的関数とした方法では，パレート最適解が得られるものの，需要家間での転売が成立し得るという意味で公平な配分になりません．
　また，目的関数の重み付けによって結果が大きく変動するという問題もあります．

## ナッシュ交渉解

　協力ゲームの分野では，ある種の公平性を満たす解としてナッシュ交渉解が知られており，各需要家の利得の対数×支払額($p_j > 0$)の和を目的関数とした以下の問題を解くことで得られます．[^1][^2]

[^1]:ナッシュ交渉解はより一般の問題で定義することが可能ですが，本記事では，資源配分問題おいてナッシュ交渉解をどう解釈できるかを紹介します．

[^2]:ゲーム理論の分野でナッシュ交渉解を扱うときには目的関数を利得の積とすることが多いのですが，後で解くときに扱いやすいように，目的関数を対数の和としています．また，$y_j = 0$ とならないようにこっそり $y_j > 0$ としています．

```math
\begin{align}
\text{max} \ & \ \sum_{j=1}^{n} p_j \log(y_j) \\
\text{s.t.} \ & \ \sum_{i=1}^m w_{ij} x_{ij} = y_j \\
& \ \sum_{j=1}^n x_{ij} \leq s_i \\
& \ x_{ij} \geq 0 \\
& \ y_j > 0
\end{align}
```

　ナッシュ交渉解の数学的な性質に関してはインターネット上で確認できる資料などでも紹介されているので，本記事ではナッシュ交渉解をどのように解釈できるかを紹介します．
 
　これまでの話でほぼネタバレしていますが，今扱っている資源配分問題において，ナッシュ交渉解は
 
:::note
パレート最適かつ，転売が成立しない均衡状態を達成する解である
:::

と解釈することができます．

　$p_j \log(y_j)$の単調性からパレート最適解が得られることは明らかなので，ナッシュ交渉解では転売が成立しないことを確認していきます．

　さきほどの具体例でナッシュ交渉解を計算すると，カレーライス1食分とラーメン0.4食分がエネルギッシュに配分され利得は 800kcal，ラーメン0.6食分とプロテイン1食分がマッスルに配分され利得は 32g という結果が得られます．


![graph3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/ca88d7f3-0cca-457a-4b47-c19ea4bd3275.png)

　利得の線形和の場合と同様に各需要家にとっての資源の単価を計算すると以下のようになります．

![image5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/9ec45241-f091-e92f-1ca0-9cf364659fb9.png)

　エネルギッシュとマッスルにとってのラーメンの単価はどちらも￥625なので，エネルギッシュとマッスルの間ではラーメンの転売が成立しません．また，カレーライスは全てエネルギッシュに配分されていますが，カレーライスの単価を最も高く評価しているのがエネルギッシュであるため，エネルギッシュから他の需要家への転売は成立しません．プロテインについても同様です．

<details><summary><b>【参考】マッスルの支払額を￥1100としたケース</b></summary>

　ちゃんと転売が成立しない解が得られています．また，利得の線形和を目的関数とした場合とは異なり，利得の変動は納得感のある量となっています．

![image6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/581d5909-b45e-7574-0eae-a3235caf801a.png)

</details>

<details><summary><b>【参考】需要家を増やしたケース</b></summary>

　ラーメンだけを欲するラーメンマンを追加します．ラーメンマンの支払額は￥500としました．需要家を増やしても，やはり転売が成立しない解が得られています．

![image7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/94e71bdb-84f0-ce9b-e1dd-6317e3f1cc16.png)

</details>

<details><summary><b>【参考】転売が成立しないことの証明</b></summary>

　 $x^\*, y^\*$ を上記の問題の最適解とし， $x^\*, y^\*$ において転売が成立しないことを背理法を使って示します．

　ある資源 $i$ と2人の需要家 $j_1, j_2$ が存在し，需要家 $j_1$ から $j_2$ への資源 $i$ の転売が成立していると仮定する．転売が成立するとは，資源 $i$ が需要家 $j_1$ に配分されており，かつ需要家 $j_1$ にとっての資源 $i$ の単価が需要家 $j_2$ にとっての資源 $i$ の単価よりも小さい
```math
x^*_{ij_1} > 0 \land \frac{p_{j_1} w_{ij_1}}{y^*_{j_1}} < \frac{p_{j_2}w_{ij_2}}{y^*_{j_2}}
```
ということである．
　ここで，任意の $t \in [0, x^*_{ij_1})$に対して
 ```math
\begin{align}
x'_{ij_1} &=& x^*_{i{j_1}} - t \\
x'_{ij_2} &=& x^*_{i{j_2}} + t \\
y'_{j_1} &=& y^*_{j_1} - w_{ij_1} t \\
y'_{j_2} &=& y^*_{j_2} + w_{ij_2} t
\end{align}
```
は実行可能解であり， $t$ についての目的関数の微分
```math
\frac{d}{d t} \sum_{j=1}^{n} p_j \log(y'_j) = -\frac{p_{ij_1} w_{ij_1}}{y^*_{j_1} - w_{ij_1} t} + \frac{p_{ij_2} w_{ij_2}}{y^*_{j_2} + w_{ij_2} t}
```
は $t=0$ のときに(単価の不等式から)正であるため， $t$ を正の範囲で十分に小さくとると
```math
\sum_{j=1}^{n} p_j \log(y'_j) > \sum_{j=1}^{n} p_j \log(y^*_j)
```
となり，これは $x^\*, y^\*$ が最適解であることと矛盾．
　よって転売は成立しない．

</details>

<details><summary><b> 【参考】供給ロジックの一般化 </b></summary>

　ところで，転売が成立しないことの証明では供給量の制約条件を使っていません．[^3]
[^3]:あらかた記事を書いたあとに気づきました．結論はわかっているからといってさぼらず，まじめに書いてみるものですね...

　今回の問題設定では各資源 $i$ の供給量を定数 $s_i$ としていますが，実は資源の供給量を決める仕組みは何でもよく，以下のような状況にも適用可能です．

* 限られた予算で資源を購入し，購入した資源を需要家に配分する
* 供給できる資源の種類の数に制限がある

 ![image8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/33f10dd7-0500-8188-323f-4ae829409295.png)

</details>

　以上のように，ナッシュ交渉解はパレート最適解であり，かつ需要家間での転売が成立しないという意味で公平な配分になっています．

# LP ソルバーでナッシュ交渉解を計算する

　ナッシュ交渉解を求めることで最適かつ公平な資源配分を実現できることがわかりましたが，目的関数が対数関数の和であり，非線形最適化問題になっています．非線形最適化のソルバーを使えるならよいのですが，諸事情で使えないこともあるので，LPソルバーを使って近似解を求める方法を紹介します．

　近似のための準備として，中間変数 $z$ を使ってナッシュ交渉解を求める問題を以下のように変形します．

```math
\begin{align}
\text{max} \ & \ \sum_{j=1}^{n} p_j z_j \\
\text{s.t.} \ & \ z_j \leq \log(y_j) \\
& \ \sum_{i=1}^m w_{ij} x_{ij} = y_j \\
& \ \sum_{j=1}^n x_{ij} \leq s_i \\
& \ x_{ij} \geq 0 \\
& \ y_j > 0 \\
& \ z_j \in \mathbb{R}
\end{align}
```

　この問題は最大化問題かつ$p_j$が正で最適解においては$z_j = \log(y_j)$となるので，もとの問題と等価であることがわかります．
　
　非線形項を含む制約条件 $z_j \leq \log(y_j)$ を，複数の線形制約

$$z_j \leq a_{jk} y_j + b_{jk}, \ k=1, \cdots $$

で近似します．

![graph_log.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635638/e8005790-995e-1bc6-5541-6b1ee8dbe466.png)

　近似解は以下の線形計画問題を解くことで得られます．

```math
\begin{align}
\text{max} \ & \ \sum_{j=1}^n p_j z_j \\
\text{s.t.} \ & \ z_j \leq a_{jk} y_j + b_{jk} \\
& \ \sum_{i=1}^m w_{ij} x_{ij} = y_j \\
& \ \sum_{j=1}^n x_{ij} \leq s_i \\
& \ x_{ij} \geq 0 \\
& \ y_j > 0 \\
& \ z_j \in \mathbb{R}
\end{align}
```

　線形制約 $z_j \leq a_{jk} y_j + b_{jk}$ の生成方法には，

* 事前に決められた数の制約を生成しておくか，逐次的に生成するか
* （生成された制約条件の実行可能領域の）頂点が $z=\log(y)$ に接するのか，辺が接するのか，どちらでもないのか

など，さまざまな方法が考えられます．

　今回は $z_j=\log(y_j)$ の接平面(2次元なので接線)を逐次的に生成する切除平面法を使いたいと思います．事前に制約条件を生成する方法と異なり，反復を繰り返すことで高い精度の近似解を得ることが可能です．

**【切除平面法のアルゴリズム】**

1. 各 j に対して初期の制約条件 2 つを生成
    * 十分小さい $y_j$ で $z_j=\log(y_j)$ に接する線形制約
        * $y_j$ が 0 に張り付くのを防ぐため
    * 十分大きい $y_j$ で $z_j=\log(y_j)$ に接する線形制約
        * $z_j$ が極端に大きい値となることを防ぐため
2. 線形計画問題を解く(最適解を $x^\*, y^\*, z^\*$ とする)
3. $z^\*$ と $\log(y^\*) $ を比較し，十分に近い値が得られていれば終了
4. 各 j に対して， $y^\*_j$ で $z_j=\log(y_j)$ に接する線形制約[^4]を生成して問題に追加
5. 2 へ戻る

[^4]:この方法が最も簡単に実装可能と思われますが，より収束の速い切除平面の作り方はあるかもしれません．

**【コード】**

<details><summary>model.py</summary>

```python
from pydantic import BaseModel
import math
import pulp


class Consumer(BaseModel, frozen=True):
    """需要家"""

    id: str


class Resource(BaseModel, frozen=True):
    """資源"""

    id: str


class Model:

    MIN_GAIN = 1e-4
    """需要家が得る価値の下限"""


    class Result(BaseModel, frozen=True):
        """結果"""

        allocation: dict[tuple[Resource, Consumer], float]
        """資源の需要家への割当量"""

        gain: dict[Consumer, float]
        """需要家の利得"""

        pricing: dict[tuple[Resource, Consumer], float]
        """資源 1 単位の需要家にとっての金額"""


    def __init__(
        self,
        supply: dict[Resource, float],
        payment: dict[Consumer, float],
        weight: dict[tuple[Resource, Consumer], float],
    ):
        """
        モデルを構築

        :param supply: 資源の供給量
        :param payment: 需要家の支払額
        :param weight: 資源・需要家ごとの，資源1単位あたりの利得
        """

        assert all(s >= 0 for s in supply.values())
        assert all(p >= 0 for p in payment.values())
        assert all(w >= 0 for w in weight.values())
        assert all(i in supply and c in payment for i, c in weight.keys())

        self.supply = supply
        """資源の供給量"""

        self.payment = payment
        """需要家の支払額"""

        self.weight = weight
        """資源・需要家ごとの，資源1単位あたりの利得"""

        self.resources = list(self.supply.keys())
        """資源の集合"""

        self.consumers = list(self.payment.keys())
        """需要家の集合"""

        self.P = pulp.LpProblem(sense=pulp.LpMaximize)
        """最適化問題(最大化)"""

        # 変数

        self.x: dict[(Resource, Consumer), pulp.LpVariable] = {}
        """配分量"""
        for (resource, consumer), weight in self.weight.items():
            if weight > 0:
                self.x[resource, consumer] = pulp.LpVariable(
                    name=f"x_{resource.id}_{consumer.id}",
                    lowBound=0,
                )

        self.y: dict[Consumer, pulp.LpVariable] = {}
        """需要家が得る利得"""
        for consumer in self.consumers:
            self.y[consumer] = pulp.LpVariable(
                name=f"y_{consumer}",
                lowBound=self.MIN_GAIN
            )

        self.z: dict[Consumer, pulp.LpVariable] = {}
        """log(y)の近似値"""
        for consumer in self.consumers:
            self.z[consumer] = pulp.LpVariable(
                name=f"z_{consumer}"
            )

        # 目的関数
        self.P += pulp.lpSum(self.payment[consumer] * self.z[consumer] for consumer in self.consumers)

        # 制約条件

        # 資源ごとの配分量の合計は，資源の供給量以下
        for resource, supply in self.supply.items():
            self.P += pulp.lpSum(
                self.x[resource, consumer] for consumer in self.consumers
                if (resource, consumer) in self.weight
            ) <= supply

        # 需要家が得る利得
        for consumer in self.consumers:
            self.P += self.y[consumer] == pulp.lpSum(
                self.weight[resource, consumer] * self.x[resource, consumer] for resource in self.resources
                if (resource, consumer) in self.weight
            )

    def solve(self, eps: float, limit_iterations: int) -> Result:
        """
        ナッシュ交渉解を計算する
        """
        print("{}\t{}\t{}\t{}\t{}".format("ITERATION", "OBJECTIVE", "UPPER", "GAP", "CUTS"), flush=True)
        # 初期カットを生成
        number_of_cuts = self.__add_initial_cuts()
        # 最大 limit_iterations 繰り返す
        for iteration in range(limit_iterations):
            # LP を解く
            self.P.solve(pulp.PULP_CBC_CMD(msg=False))
            upper = pulp.value(self.P.objective)
            # 暫定解まわりで目的関数を線形近似
            n, obj = self.__add_cuts(
                gain={consumer: pulp.value(y) for consumer, y in self.y.items()},
                eps=eps
            )
            # ログを出力
            print("{:}\t{:.7e}\t{:.7e}\t{:.7e}\t{}".format(iteration, obj, upper, upper - obj, number_of_cuts), flush=True)
            # 制約が追加されなければ(十分な精度の解が得られていれば)中断
            if n == 0:
                break
            number_of_cuts += n

        # 結果を算出
        allocation = {}
        gain = {}
        pricing = {}
        for consumer in self.consumers:
            g = pulp.value(self.y[consumer])
            gain[consumer] = g
            for item in self.resources:
                allocation[item, consumer] = pulp.value(self.x[item, consumer]) if (item, consumer) in self.x else 0
                pricing[item, consumer] = self.payment[consumer] * self.weight.get((item, consumer), 0) / g

        return self.Result(allocation=allocation, gain=gain, pricing=pricing)


    def __add_initial_cuts(self) -> int:
            # MIN_EVALUATION まわりで目的関数を線形近似
            n1, _ = self.__add_cuts(
                gain={consumer: self.MIN_GAIN for consumer in self.consumers},
                eps=0
            )
            # 需要家の利得の最大値を計算
            max_gain = {consumer: 0 for consumer in self.consumers}
            for (resource, consumer), weight in self.weight.items():
                max_gain[consumer] += weight * self.supply[resource]
            # 需要家の利得の最大値まわりで目的関数を線形近似
            n2, _ = self.__add_cuts(
                gain=max_gain,
                eps=0
            )
            return n1 + n2

    def __add_cuts(self, gain: dict[Consumer, float], eps: float) -> tuple[int, float]:
        # 追加した制約の本数
        n = 0
        # 真の目的間数値
        obj = 0
        # consumer ごとに制約を追加する
        for consumer, y0 in gain.items():
            assert y0 > 0
            # 緩和した目的関数値
            z_relaxed = pulp.value(self.z[consumer]) if pulp.value(self.z[consumer]) is not None else math.inf
            # y0 での真の目的関数値
            z0 = math.log(y0)
            assert z0 <= z_relaxed + 1e-6  # z0 は z_relaxed 以下であるはず
            # y0 での目的関数の勾配
            dz_dy = 1 / y0
            # ギャップが eps より大きければ (y0, z0) で線形近似した制約を追加
            if z0 < z_relaxed - eps:
                self.P += self.z[consumer] <= z0 + dz_dy * (self.y[consumer] - y0)
                n += 1
            obj += self.payment[consumer] * z0
        return n, obj
```

</details>

<details><summary>examples.py</summary>

```python
import random
import math

from model import Consumer, Resource, Model


def example1():
    supply = {
        Resource(id="カレー"): 1,
        Resource(id="ラーメン"): 1,
        Resource(id="プロテイン"): 1,
    }
    payment = {
        Consumer(id="エネルギッシュ"): 1000,
        Consumer(id="マッスル"): 1000,
    }
    weight = {
        (Resource(id="カレー"), Consumer(id="エネルギッシュ")): 600,
        (Resource(id="カレー"), Consumer(id="マッスル")): 10,
        (Resource(id="ラーメン"), Consumer(id="エネルギッシュ")): 500,
        (Resource(id="ラーメン"), Consumer(id="マッスル")): 20,
        (Resource(id="プロテイン"), Consumer(id="エネルギッシュ")): 100,
        (Resource(id="プロテイン"), Consumer(id="マッスル"),): 20,
    }
    model = Model(supply=supply, payment=payment, weight=weight)
    result = model.solve(1e-7, 100)
    print("GAIN")
    for consumer, gain in result.gain.items():
        print(consumer.id, gain)
    print("ALLOCATION, PRICING")
    for resource, consumer in result.allocation:
        print(resource.id, consumer.id, result.allocation[resource, consumer], result.pricing[resource, consumer])


def example2():
    supply = {
        Resource(id="カレー"): 1,
        Resource(id="ラーメン"): 1,
        Resource(id="プロテイン"): 1,
    }
    payment = {
        Consumer(id="エネルギッシュ"): 1000,
        Consumer(id="マッスル"): 1100,
    }
    weight = {
        (Resource(id="カレー"), Consumer(id="エネルギッシュ")): 600,
        (Resource(id="カレー"), Consumer(id="マッスル")): 10,
        (Resource(id="ラーメン"), Consumer(id="エネルギッシュ")): 500,
        (Resource(id="ラーメン"), Consumer(id="マッスル")): 20,
        (Resource(id="プロテイン"), Consumer(id="エネルギッシュ")): 100,
        (Resource(id="プロテイン"), Consumer(id="マッスル"),): 20,
    }
    model = Model(supply=supply, payment=payment, weight=weight)
    result = model.solve(1e-7, 100)
    print("GAIN")
    for consumer, gain in result.gain.items():
        print(consumer.id, gain)
    print("ALLOCATION, PRICING")
    for resource, consumer in result.allocation:
        print(resource.id, consumer.id, result.allocation[resource, consumer], result.pricing[resource, consumer])


def example3():
    supply = {
        Resource(id="カレー"): 1,
        Resource(id="ラーメン"): 1,
        Resource(id="プロテイン"): 1,
    }
    payment = {
        Consumer(id="エネルギッシュ"): 1000,
        Consumer(id="マッスル"): 1000,
        Consumer(id="ラーメンマン"): 500,
    }
    weight = {
        (Resource(id="カレー"), Consumer(id="エネルギッシュ")): 600,
        (Resource(id="カレー"), Consumer(id="マッスル")): 10,
        (Resource(id="ラーメン"), Consumer(id="エネルギッシュ")): 500,
        (Resource(id="ラーメン"), Consumer(id="マッスル")): 20,
        (Resource(id="ラーメン"), Consumer(id="ラーメンマン")): 1,
        (Resource(id="プロテイン"), Consumer(id="エネルギッシュ")): 100,
        (Resource(id="プロテイン"), Consumer(id="マッスル"),): 20,
    }
    model = Model(supply=supply, payment=payment, weight=weight)
    result = model.solve(1e-7, 100)
    print("GAIN")
    for consumer, gain in result.gain.items():
        print(consumer.id, gain)
    print("ALLOCATION, PRICING")
    for resource, consumer in result.allocation:
        print(resource.id, consumer.id, result.allocation[resource, consumer], result.pricing[resource, consumer])


def example_random(number_of_resources: int, number_of_consumers: int, seed: int = 42):

    rng = random.Random(seed)

    resources = [Resource(id=f"resource{i}") for i in range(number_of_resources)]
    consumers = [Consumer(id=f"consumer{j}") for j in range(number_of_consumers)]

    supply = {
        item: rng.uniform(1, 10) for item in resources
    }

    payment = {
        consumer: rng.uniform(1, 10) for consumer in consumers
    }

    weight: dict[tuple(Resource, Consumer): float] = {}
    for consumer in consumers:
        # NOTE: 全ての item に対して重みが 0 である consumer が存在するとまずいので 1 つを選んで非零にする
        weight[rng.choice(resources), consumer] = rng.uniform(1, 10)
        # それ以外は 0 または [1, 10]
        for resource in resources:
            if (resource, consumer) not in weight and rng.random() >= 0.5:
                weight[resource, consumer] = rng.uniform(1, 10)

    model = Model(supply=supply, payment=payment, weight=weight)

    result = model.solve(1e-7, 100)
    print("GAIN")
    for consumer, gain in result.gain.items():
        if gain > 1e-8:
            print(consumer.id, gain)
    print("ALLOCATION, PRICING")
    for resource, consumer in result.allocation:
        if result.allocation[resource, consumer] > 1e-8:
            print(resource.id, consumer.id, result.allocation[resource, consumer], result.pricing[resource, consumer])

    # # 転売不可能であることをチェック(重いので注意)
    # for resource in resources:
    #     for consumer1 in consumers:
    #         for consumer2 in consumers:
    #             if result.pricing.get((resource, consumer1), -math.inf) * 1.01 < result.pricing.get((resource, consumer2), math.inf):
    #                 # consumer1 よりも consumer 2 の方が resource を高く評価している ⇒ consumer1 には resource が配分されていない
    #                 assert result.allocation[resource, consumer1] <= 1e-6
    # print("SUCCESS")


example1()

example2()

example3()

# example_random(100, 100)
```

</details>

　これまでに挙げた具体例でのナッシュ交渉解はこのコードで計算しています．
　また，ランダムに生成したデータで実験したところ，「100 資源 × 100 需要家」程度の規模であれば数秒で解くことができるようです．

# まとめ

* 以下の前提で，資源を需要家に対して最適かつ公平に配分する問題を扱った
    * 資源の配分量は連続
    * 需要家の利得関数は線形
* 解を求める 2 つの方法を紹介
    * 利得の線形和を目的関数とする方法
        * パレート最適解が得られる
        * 需要家間での転売が成立するという意味で公平な配分にならない
    * ナッシュ交渉解
        * パレート最適解が得られる
        * 需要家間での転売が成立しないという意味で公平な配分が得られる
* ナッシュ交渉解をLPソルバーで求める方法を紹介

# あとがき

　何年か前に資源配分の公平性について考えていて「転売が成立しない均衡状態」を思いつき，ナッシュ交渉解の面白い解釈だと思ったので記事にしました．

　実務で使う機会がなく実用性のほどはよくわかっていないのですが，今回実験を行ったことで，ある程度までの問題規模であれば実用的な時間で解けることを確認できました．ただ，個人的にはまだよくわかっていないこともあります．

* 利得の最大化と公平性以外の指標が存在する場合にどうモデル化するのがよいのか
    * たとえば，店舗の経済性など
* 配分量が離散である場合の公平性はどうなるのか
    * 「転売が成立しない」に近い解釈を考えることが可能なのか
    * 単に変数 $x$ を整数として，実用的に許容できる結果が得られるのか

　もし今後使う機会があって何かしらの知見が得られたら，本記事の続きを書きたいと思います．