TL;DR
Python ではじめる数理最適化（第2版） の 第２章 P 17 の問題を Go と Python で解き比べました。
問題設定
ある工場では製品 $p$ と $q$ を製造しています。製品 $p$ と $q$ を製造するには原料 $m$ と $n$ が必要で、次のことがわかっています。
- 製品 $p$ を $1$ kg 製造するには、原料 $m$ が $1$ kg、原料 $n$ が $2$ kg 必要である。
- 製品 $q$ を $1$ kg 製造するには、原料 $m$ が $3$ kg、原料 $n$ が $1$ kg 必要である。
- 原料 $m$ の在庫は $30$ kg、原料 $n$ の在庫は $40$ kg である。
- 単位量あたりの利得は、製品 $p$ は $1$ 万円、製品 $q$ は $2$ 万円である。
この条件で利得を最大にするためには、製品 $p$ と $q$ をそれぞれ何 kg 製造すればよいでしょうか。
PuLP で解いてみるよ
こちらの解説は、ぜひ書籍をお買い求めいただいてお確かめくださいー
import pulp
problem = pulp.LpProblem("LP", pulp.LpMaximize)
x = pulp.LpVariable('x', cat='Continuous')
y = pulp.LpVariable('y', cat='Continuous')
problem += 1 * x + 3 * y <= 30
problem += 2 * x + 1 * y <= 40
problem += x >= 0
problem += y >= 0
problem += x + 2 * y
status = problem.solve()
print('Status:', pulp.LpStatus[status])
print('x= ', x.value(), 'y= ', y.value(), 'obj= ', problem.objective.value())
Status: Optimal
x= 18.0 y= 4.0 obj= 26.0
Go で解いてみるよ
gonum (https://gonum.org) といういろいろ数値計算ができるライブラリがあります。
今回はその中から、lp パッケージを使用して Simplex Algorithm により、解いてみます。
package main
import (
"fmt"
"log"
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/optimize/convex/lp"
)
func main() {
c := []float64{1, 2, 0, 0}
A := mat.NewDense(2, 4, []float64{1, 3, 1, 0, 2, 1, 0, 1})
b := []float64{30, 40}
// 目的関数の最大化のため、c に -1 をかける
for i, v := range c {
c[i] = -1 * v
}
z, s, err := lp.Simplex(c, A, b, 0, nil)
if err != nil {
log.Fatal(err)
}
fmt.Printf("x= %v y= %v obj= %v\n", s[0], s[1], z*-1)
}
x= 18 y= 4 obj= 26
PuLP は Maximize で解いているのに対して、Go は Minimize で解いているので、目的関数の値が異なりますが、
等しい解を得ることができました！
Simplex Algorithm について
https://pkg.go.dev のドキュメントを読むと。。。
Strang, Gilbert. "Linear Algebra and Applications." Academic, New York (1976). の第８章に説明があって、lectures 11-13 of UC Math 352 の動画へのリンクがあります。
ふむふむ
問題文から立式するよ
$$
\displaylines{
maximize \quad 1 \cdot x + 2 \cdot y \newline
s.t. \quad x \geq 0 \newline
\quad y \geq 0 \newline
\quad 1 \cdot x + 3 \cdot y \leq 30 \newline
\quad 2 \cdot x + 1 \cdot y \leq 40 \newline
}
$$
- maximize : 目的関数
- s.t. : 制約条件
グラフで描画
標準形にする
ここで、非負の値を取るスラック変数 (slack vairables) $s_1$ と $s_2$ を導入します。
$$
\displaylines{
z = 1 \cdot x + 2 \cdot y \newline
1 \cdot x + 3 \cdot y + s_1 = 30 \newline
2 \cdot x + 1 \cdot y + s_2 = 40 \newline
x, y, s_1, s_2 \geq 0
}
$$
目的関数の式変形
目的関数 $z$ を最大化するために変形します。
$$
\displaylines{
z - 1 \cdot x - 2 \cdot y = 0 \newline
1 \cdot x + 3 \cdot y + s_1 = 30 \newline
2 \cdot x + 1 \cdot y + s_2 = 40 \newline
x, y, s_1, s_2 \geq 0
}
$$
行列で表すと
$$
\displaylines{
\boldsymbol{A} = \begin{pmatrix}
1 & 3 & 1 & 0 \newline
2 & 1 & 0 & 1
\end{pmatrix} \newline
\boldsymbol{s} = \begin{pmatrix}
x \newline
y \newline
s_1 \newline
s_2 \newline
\end{pmatrix}
\newline
\boldsymbol{b} = \begin{pmatrix}
30 \newline
40
\end{pmatrix} \newline
\boldsymbol{c} = \begin{pmatrix}
1 \newline
2 \newline
0 \newline
0 \newline
\end{pmatrix}
}
$$
として
最大化:
$$
\boldsymbol{c}^T \cdot \boldsymbol{s} =
\begin{pmatrix}
1 & 2 & 0 & 0
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y \newline
s_1 \newline
s_2 \newline
\end{pmatrix}
$$
条件:
$$
\boldsymbol{A} \cdot \boldsymbol{s} = \boldsymbol{b}, \quad \boldsymbol{s} \geq 0
$$
$$
\begin{pmatrix}
1 & 3 & 1 & 0 \newline
2 & 1 & 0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y \newline
s_1 \newline
s_2 \newline
\end{pmatrix} =
\begin{pmatrix}
30 \newline
40
\end{pmatrix}
$$
のように表すことができます。
制約条件に変数導入
制約条件にいくつかの変数を導入します。
- $\boldsymbol{s_N}$ : 非基底変数 (non-basic variables) を表す行列
- $\boldsymbol{s_B}$ : 基底変数 (basic variables) を表す行列
- $\boldsymbol{N}$ : 制約条件の係数行列 $\boldsymbol{A}$ のうち非基底変数の係数の部分行列
- $\boldsymbol{B}$ : 制約条件の係数行列 $\boldsymbol{A}$ のうち基底変数の係数の部分行列
変数の内容は以下です
$$
\displaylines{
\boldsymbol{s_N} =
\begin{pmatrix}
x \newline
y \newline
\end{pmatrix}, \quad
\boldsymbol{s_B} =
\begin{pmatrix}
s_1 \newline
s_2 \newline
\end{pmatrix} \newline
\boldsymbol{N} =
\begin{pmatrix}
1 & 3 \newline
2 & 1
\end{pmatrix}, \quad
\boldsymbol{B} =
\begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
}
$$
そうすると、$\boldsymbol{A} \cdot \boldsymbol{s}$ は次のようにあらわせられます。
$$
\displaylines{
\boldsymbol{A} \cdot \boldsymbol{s} =
\begin{pmatrix}
1 & 3 & 1 & 0 \newline
2 & 1 & 0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y \newline
s_1 \newline
s_2 \newline
\end{pmatrix} \newline
= \begin{pmatrix}
1 & 3 \newline
2 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y \newline
\end{pmatrix}
+
\begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
s_1 \newline
s_2 \newline
\end{pmatrix} \newline
= \boldsymbol{N}\cdot\boldsymbol{s_N} + \boldsymbol{B}\cdot\boldsymbol{s_B}
}
$$
よって
$$
\boldsymbol{N}\cdot\boldsymbol{s_N} + \boldsymbol{B}\cdot\boldsymbol{s_B} = \boldsymbol{b}
$$
目的関数にに変数導入
- $\boldsymbol{c_N}$ : 非基底変数の係数行列
- $\boldsymbol{c_B}$ : 基底変数の係数行列
$\boldsymbol{c}$ は次のように分解できます。
$$
\displaylines{
\boldsymbol{c_N} =
\begin{pmatrix}
1 \newline
2 \newline
\end{pmatrix}, \quad
\boldsymbol{c_B} =
\begin{pmatrix}
0 \newline
0 \newline
\end{pmatrix}
}
$$
すると、$\boldsymbol{c}^T \cdot \boldsymbol{s}$ は次のようにあらわせられます。
$$
\displaylines{
\boldsymbol{c}^T \cdot \boldsymbol{s} =
\begin{pmatrix}
1 & 2 & 0 & 0
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y \newline
s_1 \newline
s_2 \newline
\end{pmatrix} \newline
= \begin{pmatrix}
1 & 2
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
y
\end{pmatrix}
+
\begin{pmatrix}
0 & 0
\end{pmatrix}
\cdot
\begin{pmatrix}
s_1 \newline
s_2
\end{pmatrix} \newline
= \boldsymbol{c_N}^T \cdot \boldsymbol{s_N} + \boldsymbol{c_B}^T \cdot \boldsymbol{s_B}
}
$$
よって
$$
z = \boldsymbol{c_N}^T \cdot \boldsymbol{s_N} + \boldsymbol{c_B}^T \cdot \boldsymbol{s_B}
$$
ところで
以下の式の両辺に左から $\boldsymbol{B}^{-1}$ をかけると
$$
\displaylines{
\boldsymbol{N}\cdot\boldsymbol{s_N} + \boldsymbol{B} \cdot \boldsymbol{s_B} = \boldsymbol{b} \newline
\boldsymbol{B}^{-1} \cdot \boldsymbol{N}\cdot\boldsymbol{s_N} + \boldsymbol{B}^{-1} \cdot \boldsymbol{B} \cdot \boldsymbol{s_B} = \boldsymbol{B}^{-1} \cdot \boldsymbol{b} \newline
\boldsymbol{B}^{-1} \cdot \boldsymbol{N}\cdot\boldsymbol{s_N} + \boldsymbol{s_B} = \boldsymbol{B}^{-1} \cdot \boldsymbol{b} \newline
\boldsymbol{s_B} = \boldsymbol{B}^{-1} \cdot \boldsymbol{b} - \boldsymbol{B}^{-1} \cdot \boldsymbol{N}\cdot\boldsymbol{s_N} \newline
}
$$
これを先程の式に代入すると、
$$
\displaylines{
z = \boldsymbol{c_N}^T \cdot \boldsymbol{s_N} + \boldsymbol{c_B}^T \cdot \boldsymbol{s_B} \newline
= \boldsymbol{c_N}^T \cdot \boldsymbol{s_N} + \boldsymbol{c_B}^T \cdot \left( \boldsymbol{B}^{-1} \cdot \boldsymbol{b} - \boldsymbol{B}^{-1} \cdot \boldsymbol{N}\cdot\boldsymbol{s_N} \right) \newline
= \boldsymbol{c_N}^T \cdot \boldsymbol{s_N} + \boldsymbol{c_B}^T \cdot \boldsymbol{B}^{-1} \cdot \boldsymbol{b} - \boldsymbol{c_B}^T \cdot\boldsymbol{B}^{-1} \cdot \boldsymbol{N} \cdot \boldsymbol{s_N} \newline
= \boldsymbol{c_B}^T \cdot \boldsymbol{B}^{-1} \cdot \boldsymbol{b} \left( \boldsymbol{c_N}^T - \boldsymbol{c_B}^T \cdot\boldsymbol{B}^{-1} \cdot \boldsymbol{N} \right) \boldsymbol{s_N} \newline
}
$$
ここで
$$
\boldsymbol{s_N} = \begin{pmatrix}
x \newline
y
\end{pmatrix}
= \begin{pmatrix}
0 \newline
0
\end{pmatrix}
$$
と固定すると、
$$
\displaylines{
\boldsymbol{s_B} = \begin{pmatrix}
s_1 \newline
s_2
\end{pmatrix}
= \boldsymbol{B}^{-1} \cdot \boldsymbol{b} \newline
= \begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
30 \newline
40
\end{pmatrix} \newline
= \begin{pmatrix}
30 \newline
40
\end{pmatrix}
}
$$
となり、$\boldsymbol{s_B}$ が求まります。
また、
- $\boldsymbol{b}^\prime = \boldsymbol{B}^{-1} \cdot \boldsymbol{b}$
- $\boldsymbol{c_N}^\prime = \boldsymbol{c_N}^T - \boldsymbol{c_B}^T \cdot\boldsymbol{B}^{-1} \cdot \boldsymbol{N}$
- $\boldsymbol{N}^\prime = \boldsymbol{B}^{-1} \cdot \boldsymbol{N}$
とおくと、
- $\boldsymbol{z} = \boldsymbol{c_B}^T \cdot \boldsymbol{b}^\prime + \boldsymbol{c_N}^\prime \cdot \boldsymbol{s_N}$
- $\boldsymbol{s_B} = \boldsymbol{b}^\prime - \boldsymbol{N}^\prime \cdot \boldsymbol{s_N}$
となります。
$\boldsymbol{c_N}^\prime$ は非基底変数の係数ベクトルであり、被約費用とよばれます。
$\boldsymbol{c_N}^\prime$ の各要素は、対応する非基底変数の値を $1$ 増加させたときに目的変数の改善料を表します。
ここで、$\boldsymbol{c_N}^\prime$ の各要素を $c_1^\prime, \quad c_2^\prime$ とし、対応する非基底変数を $x, y$ と置きます。
このとき、目的関数の係数が正となる非基底変数が存在する場合、すなわち $c_k^\prime > 0$ となる $c_1, \quad c_2, \quad x, \quad y$ が存在する場合は、
$x, \quad y$ を増加させることで、目的関数の値が改善します。
非基底変数のいずれかを選んで、選ばなかった日基底変数を $0$ に固定したまま、選んだ非基底変数を増加させます。
ここで、$\boldsymbol{a_y}$ を非基底変数 $y$ に対応する $\boldsymbol{N}^\prime$ の列ベクトルとすると
$$
\displaylines{
\boldsymbol{a_y} = \begin{pmatrix}
\end{pmatrix}
3 \newline
1
}
$$
となります。
非基底変数 $y$ の値を $ y = \theta$ まで増加させた場合の目的関数と基底変数はそれぞれ、
$$
\displaylines{
\boldsymbol{z} = \boldsymbol{c_B}^T \cdot \boldsymbol{b}^\prime + c_2^\prime \cdot \theta \newline
= (0 0) \cdot
\begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
30 \newline
40
\end{pmatrix}
+
2 \theta \newline
= 2 \theta
}
$$
$$
\displaylines{
\boldsymbol{s_B} = \boldsymbol{b}^\prime - \theta \cdot \boldsymbol{a_y}
}
$$