数理モデルを使って最適化問題を解いてみる

ちょっと仕事で見かけたLpProblemってのが面白そうだったので、自分でも試してみることにした。

# 最適化問題とはなんぞ？
このへん大学の時にやったはずなのだが、きれいさっぱりだった。

https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721

# 今回の想定ケース
以下のような条件を想定。

>工場が複数、各工場で複数の商品を製造、いくつかの顧客から製品のオーダーがくるので、各工場に対する割り当てを最適化させたい。

サンプルとして以下のデータを想定

|工場|作れるもの|
|-----------|------------------|
|工場A|製品a, 製品b, 製品c|
|工場B|製品a, 製品b, 製品d|
|工場C|製品a, 製品c, 製品d|

- 需要

|顧客|欲しいもの|
|-----------|------------------|
|顧客１|製品aを10個、製品bを10個|
|顧客２|製品aを5個、製品bを2個、製品cを10個|
|顧客３|製品aを10個、製品dを10個|
|顧客４|製品aを5個、製品bを5個、製品cを5個、製品dを5個|

## 定式化

https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721#%E5%AE%9A%E5%BC%8F%E5%8C%96

どういう風なモデルにすると解けるのかを検討する

- 目的
    - それぞれの製品をどこの工場で生産したらいいかを決めたい
- どうなると嬉しいのか？
    - 全工場の生産コストが一番小さいと嬉しい
- 制限事項
    - 各工場キャパシティ？：とりあえずは無しとしておく

上記をどの典型問題に当てはめればいいのかを考える。
顧客と工場という関係がでてきたので、最初は「経路問題」になるかと思ったが、顧客はとりあえず無視できるので、単純に割当問題でできそうな感じがする。
参考：https://qiita.com/SaitoTsutomu/items/0f6c1a4415d196e64314

※施設配置問題やナップサック問題にも似てる？とも思ったが、違うような感じ。この辺は何回かやったりすると慣れてくるんだろうなという感触。そんなにやらないとは思うけど。

今回は一般割当問題として解いてみることにする。
※というか、どれが最適なのかわからなかったので、実装しながら考えることにした。

# 実装試し：PuLP
参考
https://qiita.com/samuelladoco/items/703bf78ea66e8369c455
https://qiita.com/mzmttks/items/82ea3a51e4dbea8fbc17
https://docs.pyq.jp/python/math_opt/pulp.html

中でも以下の部分を主に参考にした。
https://qiita.com/samuelladoco/items/703bf78ea66e8369c455#%E4%BE%8B2

## とりあえず一回やってみる
``` python 
# 名称
factory_list = ["factory_A", "factory_B", "factory_C"]
product_list = ["product_a", "product_b", "product_c", "product_d"]
```
工場・製品の名称を決め、リストを作成

ここで最初に想定していなかった各工場の各製品に対する生産コストの概念がでてくるので適当に設定

|工場|製品a|製品b|製品c|製品d|
|-----------|------------------|------------------|------------------|------------------|
|工場A|10|10|10|13|
|工場B|11|11|11|15|
|工場C|12|12|12|20|

イメージとしては、工場Aは全体的に生産コストが低い。逆にCは全体的に高い（労組が強い的な）
製品dについてはどこの工場でも生産コストが高い（作るのがむつかしい的な）

``` python 
cc = [
      [10, 10, 10, 13],	# Aにおける各製品に対する生産コスト
      [11, 11, 11, 15],	# Bにおける各製品に対する生産コスト
      [12, 12, 12, 20],	# Cにおける各製品に対する生産コスト
     ]

c = {} # 空の辞書
for i in factory_list:
    for j in product_list:
        c[i,j] = cc[factory_list.index(i)][product_list.index(j)]
```
上記の表を項目名でアクセスできるようにする

>cc はリストであり、添え字が数値なので、
辞書 c を定義し、cc[0][0] は c["factory_A","product_a"] でアクセスできるようにする

``` python 
{('factory_A', 'product_a'): 10,
 ('factory_A', 'product_b'): 10,
 ('factory_A', 'product_c'): 10,
 ('factory_A', 'product_d'): 13,
 ('factory_B', 'product_a'): 11,
 ('factory_B', 'product_b'): 11,
 ('factory_B', 'product_c'): 11,
 ('factory_B', 'product_d'): 15,
 ('factory_C', 'product_a'): 12,
 ('factory_C', 'product_b'): 12,
 ('factory_C', 'product_c'): 12,
 ('factory_C', 'product_d'): 20}
```
こんな感じになる。

``` python 
# 数理最適化問題（最小化）を宣言
problem = pulp.LpProblem(name="Problem-2", sense=pulp.LpMinimize)
```
``` python 
# 連続変数を宣言：下限値は0としておく
x_suffixes = [(i,j) for i in factory_list for j in product_list]
x = pulp.LpVariable.dicts("x", x_suffixes, lowBound=0, cat = pulp.LpContinuous ) 
```
``` python 
# 目的関数を宣言
problem += pulp.lpSum(c[i,j] * x[i,j] for i in factory_list for j in product_list),"TotalCost"
```
算出すべきコストとしては、各工場の各製品に対するコスト×割り当てられる生産量となる

```python
# 需要：これが制限変数になる
product_demand_list = [30, 18, 15, 15]
```
```python
# 制約条件を宣言
for j in product_list:
    print("{} の需要 {}".format(j, product_demand_list[product_list.index(j)]))
    problem += pulp.lpSum(x[i,j] for i in factory_list) == product_demand_list[product_list.index(j)]
```
各工場での製品の製造数の合計が需要と一致すればいいのでこんな式にしている

```python
# 演算できるか確認
status = problem.solve()
print(pulp.LpStatus[status])
```
>Optimal

となったので、算出できることがわかる。

```python
# 出力
for i in factory_list:
    for j in product_list:
        print("{:} = {:},  "
              .format(x[i,j].name, x[i,j].value()), end="")
    print("")
```
こんな感じの結果となった。
>x_('factory_A',_'product_a') = 30.0,  x_('factory_A',_'product_b') = 18.0,  x_('factory_A',_'product_c') = 15.0,  x_('factory_A',_'product_d') = 15.0,  
x_('factory_B',_'product_a') = 0.0,  x_('factory_B',_'product_b') = 0.0,  x_('factory_B',_'product_c') = 0.0,  x_('factory_B',_'product_d') = 0.0,  
x_('factory_C',_'product_a') = 0.0,  x_('factory_C',_'product_b') = 0.0,  x_('factory_C',_'product_c') = 0.0,  x_('factory_C',_'product_d') = 0.0,  


全部の製品をAで作ればええやんってことになる。
たしかにそれもそうか。

## 制限事項を増やす
上記だと「工場Aの従業員がいっぱいいっぱいになってしまう」＆「工場A以外の工場がいらなくなってしまう」ので、
各工場の作業量？に対する上限と各工場の下限をつけることにする

## まずは下限
```python
factory_mini_list = [1,1,1]
# 各工場の製品の合計を0より上にする
for i in factory_list:
    print("{} の下限 {}".format(i, factory_mini_list[factory_list.index(i)]))
    problem += pulp.lpSum(x[i,j] for j in product_list) >= 1
```
※>が使えないので、「0より上」という条件にしたいところだが、「1以上」という条件にせざるを得ない

```python
Problem-2:
MINIMIZE
10*x_('factory_A',_'product_a') + 10*x_('factory_A',_'product_b') + 10*x_('factory_A',_'product_c') + 13*x_('factory_A',_'product_d') + 11*x_('factory_B',_'product_a') + 11*x_('factory_B',_'product_b') + 11*x_('factory_B',_'product_c') + 15*x_('factory_B',_'product_d') + 12*x_('factory_C',_'product_a') + 12*x_('factory_C',_'product_b') + 12*x_('factory_C',_'product_c') + 20*x_('factory_C',_'product_d') + 0
SUBJECT TO
_C1: x_('factory_A',_'product_a') + x_('factory_B',_'product_a')
 + x_('factory_C',_'product_a') = 30

_C2: x_('factory_A',_'product_b') + x_('factory_B',_'product_b')
 + x_('factory_C',_'product_b') = 18

_C3: x_('factory_A',_'product_c') + x_('factory_B',_'product_c')
 + x_('factory_C',_'product_c') = 15

_C4: x_('factory_A',_'product_d') + x_('factory_B',_'product_d')
 + x_('factory_C',_'product_d') = 15

_C5: x_('factory_A',_'product_a') + x_('factory_A',_'product_b')
 + x_('factory_A',_'product_c') + x_('factory_A',_'product_d') >= 1

_C6: x_('factory_B',_'product_a') + x_('factory_B',_'product_b')
 + x_('factory_B',_'product_c') + x_('factory_B',_'product_d') >= 1

_C7: x_('factory_C',_'product_a') + x_('factory_C',_'product_b')
 + x_('factory_C',_'product_c') + x_('factory_C',_'product_d') >= 1
```

### 結果
>x_('factory_A',_'product_a') = 28.0,  x_('factory_A',_'product_b') = 18.0,  x_('factory_A',_'product_c') = 15.0,  x_('factory_A',_'product_d') = 15.0,  
x_('factory_B',_'product_a') = 1.0,  x_('factory_B',_'product_b') = 0.0,  x_('factory_B',_'product_c') = 0.0,  x_('factory_B',_'product_d') = 0.0,  
x_('factory_C',_'product_a') = 1.0,  x_('factory_C',_'product_b') = 0.0,  x_('factory_C',_'product_c') = 0.0,  x_('factory_C',_'product_d') = 0.0,  

申し訳程度に働くA以外の工場。ちょっとだけでも働けばいいんやろというのが見えてくる程度の労働量である。

合計の生産コストも算出してみる

```python
合計のコスト
s = sum(c[i,j] * x[i,j] for i in factory_list for j in product_list)
```
>828.0

## 今度は上限
```python
# 上限
factory_max_list = [300,200,100]
for i in factory_list:
    print("{} の上限 {}".format(i, factory_max_list[factory_list.index(i)]))
    problem += pulp.lpSum(c[i,j] * x[i,j] for j in product_list) <= factory_max_list[factory_list.index(i)]
```
### 結果
>Infeasible

この上限だと収まらないらしいので、上限値を変更

```python
factory_max_list = [400,300,200]
```

### 結果2
>x_('factory_A',_'product_a') = 0.0,  x_('factory_A',_'product_b') = 5.5,  x_('factory_A',_'product_c') = 15.0,  x_('factory_A',_'product_d') = 15.0,  
x_('factory_B',_'product_a') = 14.772727,  x_('factory_B',_'product_b') = 12.5,  x_('factory_B',_'product_c') = 0.0,  x_('factory_B',_'product_d') = 0.0,  
x_('factory_C',_'product_a') = 15.227273,  x_('factory_C',_'product_b') = 0.0,  x_('factory_C',_'product_c') = 0.0,  x_('factory_C',_'product_d') = 0.0,  
合計のコスト：882.727273

それぞれがちょっとずつ仕事するようになった。

# まとめ
- とりあえず雰囲気はつかめたような気がする。
- これに顧客に対する輸送コストを追加するとどうなるんだろう？という疑問が出てくるので、もちっとだけやってみる。

- 今回作ったもの
    - https://github.com/RyoNakamae/LpProblem


# NP困難とは
いくつか読んでいると「NP困難」などの用語がでてくるが、今一つピンとこなかったのでこの辺を参考にした。
http://motojapan.hateblo.jp/entry/2017/11/15/082738
https://qiita.com/jkr_2255/items/e02362568ff1294fa0cd
http://d.hatena.ne.jp/nsasaki/20110414/1302787792
知らなくても上記を解くことはできるが知っているといいんだろうなとは思う。