はじめに
線形計画問題を解く課題が出るたびに、様々な大学のPDFにお世話になっているので、線形計画問題をシンプレックス法で解く際の手順を自分にとって分かりやすい書き方でまとめます。「○○する。」のような表現を用いるため、数学的な厳密さはありません。
対象
- シンプレックス法の数学的な説明は理解しているが、実際にどう解けば良いか分からない人
線形計画問題を解くことができれば良いだけの人は、おとなしくPuLPを使いましょう。シンプレックス法を理解するより、Pythonの環境構築の方が簡単だと思います。
問題
さっそく問題を考えます。
成分$p, q$を含む材料$1, 2$がある。材料$1$は$1 , \mathrm{kg}$あたり成分$p, q$をそれぞれ$16 , \mathrm{g}, , 8 , \mathrm{g}$含み、$1 , \mathrm{kg}$あたりの価格は、$2400 , 円$である。材料$2$は$1 , \mathrm{kg}$あたり成分$p, q$をそれぞれ$2 , \mathrm{g}, , 6 , \mathrm{g}$含み、$1 , \mathrm{kg}$あたりの価格は、$1800 , 円$である。材料$1, 2$を用いて、製品$S$を作ることを考える。ただし、製品$S$には$1 , \mathrm{kg}$あたり成分$p, q$がそれぞれ、$8 , \mathrm{g}, , 9 , \mathrm{g}$以上含まれるようにしたい。この時、製品$S$を作る際にかかるコストは最小でいくらか。
定式化
上の問題を式で表します。材料$1, 2$をそれぞれ$m_1, m_2 , \mathrm{kg}$使うことにすると、
\begin{align}
16 m_1 + 2 m_2 \geq 8 \\
8 m_1 + 6 m_2 \geq 9 \\
m_1, m_2 \geq 0
\end{align}
という条件の元、
z = 2400 m_1 + 1800 m_2
を最小化する問題だと考えることができます。
Python
シンプレックス法で解いた答えの確認も兼ねて、一度PuLPで解いてみます。
import pulp
prob = pulp.LpProblem('Cost', pulp.LpMinimize)
m_1 = pulp.LpVariable('m_1', 0, 1000, 'Continuous')
m_2 = pulp.LpVariable('m_2', 0, 1000, 'Continuous')
prob += 2400 * m_1 + 1800 * m_2
prob += 16 * m_1 + 2 * m_2 >= 8
prob += 8 * m_1 + 6 * m_2 >= 9
prob.solve()
print(f"Material1: {pulp.value(m_1)}")
print(f"Material2: {pulp.value(m_2)}")
print(f"Minimum : {2400 * pulp.value(m_1) + 1800 * pulp.value(m_2)}")
$python solve.py
Material1: 0.375
Material2: 1.0
Minimum : 2700.0
$m_1 = 0.375 , \mathrm{kg}, m_2 = 1.0 , \mathrm{kg}$で最小値$z = 2700 , 円$をとることが分かりました。
シンプレックス法
では、本題であるシンプレックス法の手順を示します。
Step 1. 標準形にする
もし最大化が目的であれば、目的関数に$-1$をかけ、最小化を目的とします。また、非負のスラック変数を導入し、不等式を等式にします。今回は、すでに目的が最小化なので、目的関数はそのままにして、スラック変数$m_3, m_4$を導入すると、
\begin{align}
z = 2400 m_1 + 1800 m_2 \\
16 m_1 + 2 m_2 - m_3 = 8 \\
8 m_1 + 6 m_2 - m_4 = 9 \\
m_1, m_2, m_3, m_4 \geq 0
\end{align}
のように表されます。行列で表すと、
\begin{align}
A &=
\left(
\begin{matrix}
16 & 2 & -1 & 0 \\
8 & 6 & 0 & -1
\end{matrix}
\right) \\
{\bf m} &=
\left(
\begin{matrix}
m_1, m_2, m_3, m_4
\end{matrix}
\right)^T \\
{\bf b} &=
\left(
\begin{matrix}
8 \\
9
\end{matrix}
\right) \\
{\bf c} &=
\left(
\begin{matrix}
2400 & 1800 & 0 & 0
\end{matrix}
\right)^T
\end{align}
として、
z = {\bf c}^T {\bf m} \\
A {\bf m} = {\bf b}
のように表すことができます。
Step 2. 実行可能解を求める
気合で実行可能解を求めます。もし、実行可能解をうまく求めることができない場合は、人為変数を用いる方法がありますが、今回は説明しません。今回の場合、式の数が2つなので、基底変数が2個になることに注意すると、
{\bf m} = (0, 4, 0, 15)
が実行可能解であることが分かります。この時、条件式は、
\left(
\begin{matrix}
8 & 1 & -1/2 & 0 \\
40 & 0 & -3 & 1
\end{matrix}
\right)
{\bf m} =
\left(
\begin{matrix}
4 \\
15
\end{matrix}
\right) \\
\left(
\begin{matrix}
-1200 & 0 & 900 & 0
\end{matrix}
\right)
{\bf m}
= z - 7200
と表されます。この時、${\bf r} = \left( \begin{matrix} -1200 & 0 & 900 & 0 \end{matrix} \right)$において、$m_1$の係数が負であるため、まだ最適化の余地があります。
Step 3. 基底変数を入れ替える
そこで、新たに$m_1$を基底変数に加えます。この時、前のstepで基底変数とした$m_2, m_4$のどちらかを基底変数から外す必要があリます。これは、条件式のうち、$m_1$以外を無視したものを考え,
8 m_1 = 4 \\
40 m_1 = 15
を解いて、小さくなる方に対応する基底変数を外せば良いです。この場合、
m_1 = 0.5 \\
m_1 = 0.375
となるので、下の式に対応する$m_4$を基底変数から外せば良いことが分かります。この時条件式は、
\left(
\begin{matrix}
0 & 1 & 1/10 & -1/5 \\
1 & 0 & -3/40 & 1/40
\end{matrix}
\right)
{\bf m} =
\left(
\begin{matrix}
1 \\
0.375
\end{matrix}
\right) \\
\left(
\begin{matrix}
0 & 0 & 0 & 30
\end{matrix}
\right)
{\bf m}
= z - 2700
となります。${\bf r} = \left( \begin{matrix} 0 & 0 & 0 & 30 \end{matrix} \right)$において、全ての係数が非負であるため、最適であることが分かります。もし、まだ負の係数がある場合は、もう一度Step 3を行う。今回はStep 4に進みます。
Step4. 解答する
\left(
\begin{matrix}
0 & 1 & 1/10 & -1/5 \\
1 & 0 & -3/40 & 1/40
\end{matrix}
\right)
{\bf m} =
\left(
\begin{matrix}
1 \\
0.375
\end{matrix}
\right) \\
\left(
\begin{matrix}
0 & 0 & 0 & 30
\end{matrix}
\right)
{\bf m}
= z - 2700
これらの式から、
m_1 = 0.375 \\
m_2 = 1.0 \\
のとき、
z = 2700
となることが分かります。これはPuLPを用いて求めた結果と一致します。以上が、シンプレックス法による線形計画問題を解く手順となります。