【数理最適化問題】PuLPを使った線形計画法

数理最適化問題をこれから解く。今回はその中の線形計画法を使って解いていく。
線形計画法を解くにはPuLPというライブラリがメジャーのためこれを使っていく。
数理最適化やPuLPが知りたい人はこちらを参照
PuLP による線型計画問題の解き方ことはじめ

サンプルの問題はこちらのサイトから引用している
NTTデータ数理システムNumerical Optimizer Sample問題集

配合問題

例題は次のページを参照
2.1 配合問題

鉛、亜鉛、スズの配合率が異なる9種類の合金からコストが最小で30%:30%:40%の配合になるようにするにはどうすればいいのかという問題。

配合問題

# 関数の最小化を目的として設定する
problem1 = pulp.LpProblem("Problem-1", pulp.LpMinimize)

# 変数の設定 (変数名、最小値、最大値、型)
# 9種類の合金の割合
x1 = pulp.LpVariable('X1', 0, 1, 'Continuous')
x2 = pulp.LpVariable('X2', 0, 1, 'Continuous')
x3 = pulp.LpVariable('X3', 0, 1, 'Continuous')
x4 = pulp.LpVariable('X4', 0, 1, 'Continuous')
x5 = pulp.LpVariable('X5', 0, 1, 'Continuous')
x6 = pulp.LpVariable('X6', 0, 1, 'Continuous')
x7 = pulp.LpVariable('X7', 0, 1, 'Continuous')
x8 = pulp.LpVariable('X8', 0, 1, 'Continuous')
x9 = pulp.LpVariable('X9', 0, 1, 'Continuous')

# 目的関数(最小にしたいコストを定義)
problem1 += 7.3*x1 + 6.9*x2 + 7.3*x3 + 7.5*x4 + 7.6*x5 + 6.0*x6 + 5.8*x7 + 4.3*x8 + 4.1*x9

# 制約の設定
# 合計で1(100％)とする
problem1 += x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 == 1

# 鉛の混合比
problem1 += 20*x1 + 50*x2 + 30*x3 +30*x4 + 30*x5 + 60*x6 + 40*x7 + 10*x8 + 10*x9 == 30
# 亜鉛の混合比
problem1 += 30*x1 + 40*x2 + 20*x3 +40*x4 + 30*x5 + 30*x6 + 50*x7 + 30*x8 + 10*x9 == 30
# スズの混合比
problem1 += 50*x1 + 10*x2 + 50*x3 +30*x4 + 40*x5 + 10*x6 + 10*x7 + 60*x8 + 80*x9 == 40

# 問題定義の確認
print(problem1)

# 実行
result = problem1.solve()

# 結果の確認
print("X1:" ,pulp.value(x1))
print("X2:" ,pulp.value(x2))
print("X3:" ,pulp.value(x3))
print("X4:" ,pulp.value(x4))
print("X5:" ,pulp.value(x5))
print("X6:" ,pulp.value(x6))
print("X7:" ,pulp.value(x7))
print("X8:" ,pulp.value(x8))
print("X9:" ,pulp.value(x9))
print("Cost:" ,pulp.value(problem1.objective))

定義した問題print関数で確認した際の出力。

配合問題出力1

# 問題定義の確認
Problem-1:
MINIMIZE
7.3*X1 + 6.9*X2 + 7.3*X3 + 7.5*X4 + 7.6*X5 + 6.0*X6 + 5.8*X7 + 4.3*X8 + 4.1*X9 + 0.0
SUBJECT TO
_C1: X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 = 1

_C2: 20 X1 + 50 X2 + 30 X3 + 30 X4 + 30 X5 + 60 X6 + 40 X7 + 10 X8 + 10 X9
 = 30

_C3: 30 X1 + 40 X2 + 20 X3 + 40 X4 + 30 X5 + 30 X6 + 50 X7 + 30 X8 + 10 X9
 = 30

_C4: 50 X1 + 10 X2 + 50 X3 + 30 X4 + 40 X5 + 10 X6 + 10 X7 + 60 X8 + 80 X9
 = 40

VARIABLES
X1 <= 1 Continuous
X2 <= 1 Continuous
X3 <= 1 Continuous
X4 <= 1 Continuous
X5 <= 1 Continuous
X6 <= 1 Continuous
X7 <= 1 Continuous
X8 <= 1 Continuous
X9 <= 1 Continuous

配合問題出力2

# 結果の確認
X1: 0.0
X2: 0.0
X3: 0.0
X4: 0.0
X5: 0.0
X6: 0.4
X7: 0.0
X8: 0.6
X9: 0.0
Cost: 4.98

結果を見ると6番の合金40%、8番の合金60%でコストが最小になり、4.98である。

輸送問題

例題は次のページを参照
2.2 輸送問題

2つの工場(1、2)から3つの店舗(a、b、c)へ荷物を輸送する。それぞれ供給可能量と需要量が決まっており、輸送コストを最小にする問題。

輸送問題

# 関数の最小化を目的として設定する
problem2 = pulp.LpProblem("Problem-2", pulp.LpMinimize)

# 変数の設定 (変数名、最小値、最大値、型)
# 工場から店舗への配送量を変数とし、供給量の最大を最大値とした
x1a = pulp.LpVariable('X1a', 0, 200, 'Continuous')
x1b = pulp.LpVariable('X1b', 0, 200, 'Continuous')
x1c = pulp.LpVariable('X1c', 0, 200, 'Continuous')
x2a = pulp.LpVariable('X2a', 0, 200, 'Continuous')
x2b = pulp.LpVariable('X2b', 0, 200, 'Continuous')
x2c = pulp.LpVariable('X2c', 0, 200, 'Continuous')

# 目的関数(最小にしたいコストを定義)
problem2 += 3.4*x1a + 2.2*x1b + 2.9*x1c + 3.4*x2a + 2.4*x2b + 2.5*x2c

# 制約の設定
# 各工場の供給可能量
problem2 += x1a + x1b + x1c <= 250
problem2 += x2a + x2b + x2c <= 450
# 各店舗の需要量
problem2 += x1a + x2a == 200
problem2 += x1b + x2b == 200
problem2 += x1c + x2c == 200

# 問題定義の確認
print(problem2)

# 実行
result = problem2.solve()

# 結果の確認
print("X1a:" ,pulp.value(x1a))
print("X1b:" ,pulp.value(x1b))
print("X1c:" ,pulp.value(x1c))
print("X2a:" ,pulp.value(x2a))
print("X2b:" ,pulp.value(x2b))
print("X2c:" ,pulp.value(x2c))
print("Cost:" ,pulp.value(problem2.objective))

定義した問題print関数で確認した際の出力。

輸送問題出力1

# 問題定義の確認
Problem-2:
MINIMIZE
3.4*X1a + 2.2*X1b + 2.9*X1c + 3.4*X2a + 2.4*X2b + 2.5*X2c + 0.0
SUBJECT TO
_C1: X1a + X1b + X1c <= 250

_C2: X2a + X2b + X2c <= 450

_C3: X1a + X2a = 200

_C4: X1b + X2b = 200

_C5: X1c + X2c = 200

VARIABLES
X1a <= 200 Continuous
X1b <= 200 Continuous
X1c <= 200 Continuous
X2a <= 200 Continuous
X2b <= 200 Continuous
X2c <= 200 Continuous

輸送問題出力2

# 結果の確認
X1a: -0.0
X1b: 200.0
X1c: -0.0
X2a: 200.0
X2b: 0.0
X2c: 200.0
Cost: 1620.0

工場1からは店舗bへ200
工場2からは店舗aへ200、店舗cへ200
最小コストは1620という結果がでた。
例題の答えを見るとコストは同じだが輸送量が異なる。
しかし目的(コストの最小化)は達成できているのでこの答えで問題はないと思う。

多期間計画問題

例題は次のページを参照
2.3 多期間計画問題

2種類の原料A，Bを加工して2種類の製品I, IIを生産している工場がある。向こう3ヶ月の生産計画を立てる問題。
各製品を1単位生産するために必要な原料の使用量，各製品の生産／在庫コスト，各製品の月ごとの出荷量，各原料の月ごとの利用可能量が与えられている

多期間計画問題

# 関数の最小化を目的として設定する
problem3 = pulp.LpProblem("Problem-3", pulp.LpMinimize)

# 変数の設定 (変数名、最小値、最大値、型)
# 生産(Prodution)I
Pi1 = pulp.LpVariable('Prodution I_1', 0, 170, 'Integer')
Pi2 = pulp.LpVariable('Prodution I_2', 0, 170, 'Integer')
Pi3 = pulp.LpVariable('Prodution I_3', 0, 170, 'Integer')
# 生産(Prodution)II
Pii1 = pulp.LpVariable('Prodution II_1', 0, 160, 'Integer')
Pii2 = pulp.LpVariable('Prodution II_2', 0, 160, 'Integer')
Pii3 = pulp.LpVariable('Prodution II_3', 0, 160, 'Integer')

# 在庫(Stock)I
Si1 = pulp.LpVariable('Stock I_1', 0, 170, 'Integer')
Si2 = pulp.LpVariable('Stock I_2', 0, 170, 'Integer')
# 在庫(Stock)II
Sii1 = pulp.LpVariable('Stock II_1', 0, 160, 'Integer')
Sii2 = pulp.LpVariable('Stock II_2', 0, 160, 'Integer')

# 目的関数(最小にしたいコストを定義)
problem3 += (75*Pi1 + 50*Pii1 + 8*Si1 + 7*Sii1) + (75*Pi2 + 50*Pii2 + 8*Si2 + 7*Sii2) + (75*Pi2 + 50*Pii2)

# 制約の設定
# 材料による生産数の制約
problem3 += 2*Pi1 + 7*Pii1 <= 920
problem3 += 5*Pi1 + 3*Pii1 <= 790
problem3 += 2*Pi2 + 7*Pii2 <= 750
problem3 += 5*Pi2 + 3*Pii2 <= 600
problem3 += 2*Pi3 + 7*Pii3 <= 500
problem3 += 5*Pi3 + 3*Pii3 <= 480
# 出荷数と在庫数の制約
problem3 += Pi1 - Si1 == 30
problem3 += Pii1 - Sii1 == 20
problem3 += Pi2 + Si1 - Si2 == 60
problem3 += Pii2 + Sii1 - Sii2 == 50
problem3 += Pi3 + Si2 == 80
problem3 += Pii3 + Sii2 == 90

# 問題定義の確認
print(problem3)

# 実行
result = problem3.solve()

# 結果の確認
print("---1ヶ月目---")
print("Prodution I_1:" ,pulp.value(Pi1))
print("Prodution II_1:" ,pulp.value(Pii1))
print("Stock I_1:" ,pulp.value(Si1))
print("Stock II_1:" ,pulp.value(Sii1))

print("---2ヶ月目---")
print("Prodution I_2:" ,pulp.value(Pi2))
print("Prodution II_2:" ,pulp.value(Pii2))
print("Stock I_2:" ,pulp.value(Si2))
print("Stock II_2:" ,pulp.value(Sii2))

print("---3ヶ月目---")
print("Prodution I_3:" ,pulp.value(Pi3))
print("Prodution II_3:" ,pulp.value(Pii3))
print("Cost:" ,pulp.value(problem3.objective))

定義した問題print関数で確認した際の出力。

多期間計画問題出力1

# 問題定義の確認
Problem-3:
MINIMIZE
50*Prodution_II_1 + 100*Prodution_II_2 + 75*Prodution_I_1 + 150*Prodution_I_2 + 7*Stock_II_1 + 7*Stock_II_2 + 8*Stock_I_1 + 8*Stock_I_2 + 0
SUBJECT TO
_C1: 7 Prodution_II_1 + 2 Prodution_I_1 <= 920

_C2: 3 Prodution_II_1 + 5 Prodution_I_1 <= 790

_C3: 7 Prodution_II_2 + 2 Prodution_I_2 <= 750

_C4: 3 Prodution_II_2 + 5 Prodution_I_2 <= 600

_C5: 7 Prodution_II_3 + 2 Prodution_I_3 <= 500

_C6: 3 Prodution_II_3 + 5 Prodution_I_3 <= 480

_C7: Prodution_I_1 - Stock_I_1 = 30

_C8: Prodution_II_1 - Stock_II_1 = 20

_C9: Prodution_I_2 + Stock_I_1 - Stock_I_2 = 60

_C10: Prodution_II_2 + Stock_II_1 - Stock_II_2 = 50

_C11: Prodution_I_3 + Stock_I_2 = 80

_C12: Prodution_II_3 + Stock_II_2 = 90

VARIABLES
0 <= Prodution_II_1 <= 160 Integer
0 <= Prodution_II_2 <= 160 Integer
0 <= Prodution_II_3 <= 160 Integer
0 <= Prodution_I_1 <= 170 Integer
0 <= Prodution_I_2 <= 170 Integer
0 <= Prodution_I_3 <= 170 Integer
0 <= Stock_II_1 <= 160 Integer
0 <= Stock_II_2 <= 160 Integer
0 <= Stock_I_1 <= 170 Integer
0 <= Stock_I_2 <= 170 Integer

多期間計画問題出力2

# 結果の確認
---1ヶ月目---
Prodution I_1: 98.0
Prodution II_1: 100.0
Stock I_1: 68.0
Stock II_1: 80.0
---2ヶ月目---
Prodution I_2: 8.0
Prodution II_2: 7.0
Stock I_2: 16.0
Stock II_2: 37.0
---3ヶ月目---
Prodution I_3: 64.0
Prodution II_3: 53.0
Cost: 15741.0

このような生産計画を立てればよいと結果が出力された。
例題の解答よりコストが小さくなった。例題のツールより良い結果を得てしまったかもしれない。