線形計画問題入門

1. 最適化問題とは

与えられた制約条件の下で目的関数を最小化(または最大化)する問題を総称して最適化問題と呼ぶ。目的関数は利益やコストなどを表す関数(評価関数とも呼ぶ)であり、制約条件は目的関数を最大化(最小化)する際に満たさなければならない条件である。どちらも数式で表現される。

最適化問題は、目的関数と制約条件の形によって分類される。

• 線形計画問題: 目的関数と制約条件がすべて線形(1次)式で表される。
• 非線形計画問題: 目的関数や制約条件に非線形(2次以上)の式を含む。
• 整数計画問題: 変数が整数のみを取る。

本資料では線形計画問題(linear programming problem)を扱う。整数計画問題は「整数計画問題入門」で、その応用は「ナップサック問題入門」「ビンパッキング問題入門」で扱う。

学習目標

1. 現実の問題(生産計画)から目的関数と制約条件を組み立てられる。
2. 2変数の線形計画問題を図で解き、最適解が実行可能領域の頂点にあることを説明できる。
3. PuLP で線形計画問題を解いて結果を読み取れる。

2. 生産計画問題

次のような問題を考えてみる。ある会社は2種類の製品 $P_1, P_2$ を生産している。

• 製品 $P_1, P_2$ は3種類の材料 $M_1, M_2, M_3$ からできている。
• 製品 $P_1, P_2$ を1 [kg] 生産したときの利益はそれぞれ3万円、5万円である。
• 製品を生産するのに必要な材料の量は次の表のとおりである。
• ただし材料 $M_1, M_2, M_3$ の在庫はそれぞれ140、100、120 [kg] である。

原料 \ 製品	$P_1$	$P_2$	在庫
$M_1$ [kg]	1	7	140
$M_2$ [kg]	2	4	100
$M_3$ [kg]	3	2	120
利益 [万円]	3	5	-

このとき、利益を最大化するには製品 $P_1, P_2$ をどれだけ生産すべきだろうか。

2.1 定式化

この問題の目的は「製品 $P_1, P_2$ を生産して得られる利益の最大化」であり、これが目的関数になる。目的を達成するための条件は「材料の使用量が在庫を超えないこと」であり、これが制約条件になる。目的関数と制約条件を数式で表現することを定式化と呼ぶ。

まず決定変数を定める。$x_1, x_2$ を製品 $P_1, P_2$ の生産量 [kg] とする。目的関数(総利益)は、1 kg あたりの利益に生産量を掛けて足し合わせた $3x_1 + 5x_2$ である。制約条件は材料ごとに1本ずつ立てる。例えば材料 $M_1$ は、$P_1$ に1 kg、$P_2$ に7 kg 使われるから、使用量は $x_1 + 7x_2$ であり、これが在庫140以下でなければならない。生産量は負にならないので非負条件も加える。まとめると次のようになる。

$$ \begin{align} \max \ \ & z = 3x_1 + 5x_2 \\ \text{s.t. } \ \ & x_1 + 7x_2 \leq 140, \\ & 2x_1 + 4x_2 \leq 100, \\ & 3x_1 + 2x_2 \leq 120, \\ & x_1 \geq 0, \ x_2 \geq 0. \end{align} $$

s.t. はsubject to の略で、制約条件と同義である。目的関数と制約条件がすべて1次式なので、これは線形計画問題である。

3. 図で解く

変数が2つなので、この問題は平面上に描いて解ける。各制約条件の不等号を等号に直すと直線になり、不等式はその直線の片側の領域を表す。すべての制約(非負条件を含む)を満たす領域、すなわち解の候補となる点の集合を実行可能領域と呼ぶ。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.style.use('default')
plt.rcParams['mathtext.fontset'] = 'cm'

x1 = np.arange(0, 200, 0.1)
x2_1 = -(1 / 7) * x1 + 20     # x1 + 7x2 = 140
x2_2 = -(1 / 2) * x1 + 25     # 2x1 + 4x2 = 100
x2_3 = -(3 / 2) * x1 + 60     # 3x1 + 2x2 = 120

y1 = np.zeros_like(x1)
y2 = np.minimum(np.minimum(x2_1, x2_2), x2_3)

plt.plot(x1, x2_1, zorder=0, label=r'$x_1+7x_2=140$')
plt.plot(x1, x2_2, zorder=0, label=r'$2x_1+4x_2=100$')
plt.plot(x1, x2_3, zorder=0, label=r'$3x_1+2x_2=120$')
plt.fill_between(x1, y1, y2, where=y1 < y2, facecolor='yellow', alpha=0.5)

plt.xlim(0, 60)
plt.ylim(0, 30)
plt.xlabel(r'$x_1$', fontsize=13)
plt.ylabel(r'$x_2$', fontsize=13)
plt.legend(fontsize=11)
plt.tight_layout()
plt.show()

黄色の領域が実行可能領域である。直線で囲まれた凸な多角形になっている。

次に目的関数 $z = 3x_1 + 5x_2$ を考える。$z$ の値を1つ固定すると、この式は1本の直線(等高線)になる。$z$ を大きくしていくと等高線は右上へ平行移動していくから、「実行可能領域と交わったまま、等高線をどこまで右上に動かせるか」を考えれば最大値が求まる。

plt.plot(x1, x2_1, zorder=0, label=r'$x_1+7x_2=140$')
plt.plot(x1, x2_2, zorder=0, label=r'$2x_1+4x_2=100$')
plt.plot(x1, x2_3, zorder=0, label=r'$3x_1+2x_2=120$')
plt.fill_between(x1, y1, y2, where=y1 < y2, facecolor='yellow', alpha=0.5)

# 目的関数の等高線 z = 60, 100, 142.5
for i, z in enumerate([60, 100, 142.5]):
    lbl = r'$z=3x_1+5x_2$' if i == 0 else None
    plt.plot(x1, -(3 / 5) * x1 + z / 5, color='red', linestyle='dotted', label=lbl)

# 実行可能領域の頂点
verts = [(0, 0), (40, 0), (35, 7.5), (14, 18), (0, 20)]
plt.scatter([v[0] for v in verts], [v[1] for v in verts], color='red', zorder=3)
for name, (vx, vy) in zip(['O', 'A', 'B', 'C', 'D'], verts):
    plt.text(vx + 1, vy + 0.7, name, fontsize=12)

plt.xlim(0, 60)
plt.ylim(0, 30)
plt.xlabel(r'$x_1$', fontsize=13)
plt.ylabel(r'$x_2$', fontsize=13)
plt.legend(fontsize=11)
plt.tight_layout()
plt.show()

等高線を右上へ動かしていくと、実行可能領域から離れる直前に残るのは領域の「角」である。線形計画問題では、最適解のうち少なくとも1つは実行可能領域の頂点に存在する。したがって頂点をすべて調べれば最適解が見つかる。頂点は2本以上の直線が交わる点だから、連立方程式を解けば求まる。

# 頂点の列挙と目的関数値の計算
z = lambda x1v, x2v: 3 * x1v + 5 * x2v
print('頂点 | (x1, x2) | z')
for name, (vx, vy) in zip(['O', 'A', 'B', 'C', 'D'], verts):
    print(f'{name} | ({vx}, {vy}) | {z(vx, vy)}')

頂点 | (x1, x2) | z
O | (0, 0) | 0
A | (40, 0) | 120
B | (35, 7.5) | 142.5
C | (14, 18) | 132
D | (0, 20) | 100

頂点Bの $(x_1, x_2) = (35, 7.5)$ で最大値 $z = 142.5$ となる。つまり $P_1$ を35 kg、$P_2$ を7.5 kg 作れば最大利益142.5万円が得られる。なお頂点Bは直線 $2x_1 + 4x_2 = 100$ と $3x_1 + 2x_2 = 120$ の交点であり、連立方程式を解いて確かめられる。

4. PuLP で解く

図で解けるのは変数が2つのときだけである。実務では変数が100個、制約が50本といった問題を扱うため、ソルバー(最適化問題を解くソフトウェア)を使う。PythonからはPuLP ライブラリを通じてソルバーを呼び出せる。

なおGoogle Colab で実行する際には、次のセルの2行目のコメントを外して実行すること。

# Google Colab の場合は2行目のコメントを外して実行: ライブラリインストール(1回だけ)
# !pip install pulp matplotlib

import pulp

prob = pulp.LpProblem('production_planning', pulp.LpMaximize)
x1v = pulp.LpVariable('x1', lowBound=0)  # lowBound=0 が非負条件
x2v = pulp.LpVariable('x2', lowBound=0)

prob += 3 * x1v + 5 * x2v          # 目的関数
prob += x1v + 7 * x2v <= 140       # 材料M1
prob += 2 * x1v + 4 * x2v <= 100   # 材料M2
prob += 3 * x1v + 2 * x2v <= 120   # 材料M3

prob.solve(pulp.PULP_CBC_CMD(msg=0))
print('求解結果:', pulp.LpStatus[prob.status])
print(f'最適解: x1 = {x1v.value()}, x2 = {x2v.value()}')
print('最適値: z =', pulp.value(prob.objective))

求解結果: Optimal
最適解: x1 = 35.0, x2 = 7.5
最適値: z = 142.5

図で求めた頂点Bと同じ $(35, 7.5)$、$z = 142.5$ が得られた。

ソルバーの内部では単体法(simplex method)と呼ばれるアルゴリズムが使われている。単体法は「最適解は頂点にある」という線形計画問題の性質を利用し、ある頂点から出発して目的関数が最も増える隣の頂点へ移動することを繰り返す。すべての頂点を調べる必要がなく、変数が増えても機械的に適用できるため、大規模な問題でも現実的な時間で解ける。単体法の詳しい手順は授業で扱う。

5. 課題

演習問題1

あるスーパーはお惣菜 $D_1, D_2$ を材料 $M_1, M_2, M_3$ を使って生産している。必要な材料の量と1 kg あたりの利益、各材料の在庫量は次の表のとおりである。利益を最大化するにはお惣菜 $D_1, D_2$ をどれだけ生産すればよいか。

材料 \ お惣菜	$D_1$	$D_2$	在庫
$M_1$ [kg]	1	2	14
$M_2$ [kg]	1	1	8
$M_3$ [kg]	3	1	18
利益 [万円/kg]	2	3	-

1. この問題を定式化せよ。変数 $x_1, x_2$ をそれぞれお惣菜 $D_1, D_2$ の生産量とし、生産量は必ず0以上であるとする。
2. 定式化した問題の制約条件を等式になおし、2次元グラフ上に図示せよ。
3. 図示したグラフ上で、この問題の実行可能領域を斜線で塗りつぶせ。
4. 実行可能領域の各頂点について、変数の値と目的関数の値を求めよ。
5. 最適解とそのときの目的関数の値を示せ。

演習問題2

お惣菜の材料が $M_4, M_5, M_6$ に変更された。それに伴い、材料の使用量や在庫、お惣菜の利益が次の表のように変わった。演習問題1と同じ手順(定式化、図示、頂点の列挙、最適解)で、利益を最大化する生産量を求めよ。

材料 \ お惣菜	$D_1$	$D_2$	在庫
$M_4$ [kg]	3	1	9
$M_5$ [kg]	12	7	42
$M_6$ [kg]	4/3	2	8
利益 [万円/kg]	3	2	-

問題文の原本は同フォルダの LP-prob.docx にある。解答は各自でWord ファイルにまとめて提出すること。

6. 注意事項

• 定式化に慣れないうちは、表の1行が制約1本に対応することを意識するとよい。「何が変数か」「何を最大化(最小化)するか」「守るべき条件は何か」の3つを、数式にする前に日本語で書き出すこと。
• 手で解いた結果とPuLP の結果が一致することを必ず確かめること。一致しなければ定式化かプログラムのどちらかが間違っている。
• 理解できない場合には指導教員まで相談すること。