Edited at

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

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


最適化問題とはなんぞ?

このへん大学の時にやったはずなのだが、きれいさっぱりだった。

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


今回の想定ケース

以下のような条件を想定。


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


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

工場
作れるもの

工場A
製品a, 製品b, 製品c

工場B
製品a, 製品b, 製品d

工場C
製品a, 製品c, 製品d


  • 需要

顧客
欲しいもの

顧客1
製品aを10個、製品bを10個

顧客2
製品aを5個、製品bを2個、製品cを10個

顧客3
製品aを10個、製品dを10個

顧客4
製品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


とりあえず一回やってみる

# 名称

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についてはどこの工場でも生産コストが高い(作るのがむつかしい的な)

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"] でアクセスできるようにする


{('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}

こんな感じになる。

# 数理最適化問題(最小化)を宣言

problem = pulp.LpProblem(name="Problem-2", sense=pulp.LpMinimize)

# 連続変数を宣言:下限値は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 )

# 目的関数を宣言

problem += pulp.lpSum(c[i,j] * x[i,j] for i in factory_list for j in product_list),"TotalCost"

算出すべきコストとしては、各工場の各製品に対するコスト×割り当てられる生産量となる

# 需要:これが制限変数になる

product_demand_list = [30, 18, 15, 15]

# 制約条件を宣言

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)]

各工場での製品の製造数の合計が需要と一致すればいいのでこんな式にしている

# 演算できるか確認

status = problem.solve()
print(pulp.LpStatus[status])


Optimal


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

# 出力

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以外の工場がいらなくなってしまう」ので、

各工場の作業量?に対する上限と各工場の下限をつけることにする


まずは下限

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以上」という条件にせざるを得ない

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以外の工場。ちょっとだけでも働けばいいんやろというのが見えてくる程度の労働量である。

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

合計のコスト

s = sum(c[i,j] * x[i,j] for i in factory_list for j in product_list)


828.0



今度は上限

# 上限

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


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

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


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


まとめ


  • とりあえず雰囲気はつかめたような気がする。

  • これに顧客に対する輸送コストを追加するとどうなるんだろう?という疑問が出てくるので、もちっとだけやってみる。



  • 今回作ったもの




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

知らなくても上記を解くことはできるが知っているといいんだろうなとは思う。