概要
線形計画問題の図解法を解説するための資料を matplotlib （Python）で作成するためのメモです。講義資料や模範解答の作成などに・・・。
問題例
次のような線形計画問題（LP）を考えます。$\mathbb{R}$ は実数集合の意味で使っています。
線形計画問題(p.1)
\begin{align}
\mathrm{Minimize}\quad & 100x_1 + 300x_2 & \quad (0) \\
\mathrm{subject}\;\mathrm{to}\quad & 4x_1 + x_2 \geq -22 & \quad (1)\\
& x_1 + 2x_2 \geq -13 & \quad (2)\\
& x_1 + 5x_2 \geq -29 & \quad (3)\\
& x_2 \leq -1 & \quad (4)\\
& x_1,\;x_2 \in \mathbb{R }& \quad (5)
\end{align}
求解
線形計画問題の求解自体は scipy.optimize や、その他のライブラリで簡単にできます。scipy.optimize で求解するには、まず、次のような形に問題を整理します（不等号の向きに注意）。
線形計画問題(p.1)
\begin{align}
\mathrm{Minimize}\quad & \boldsymbol{c}^T \boldsymbol{x}\\
\mathrm{subject}\;\mathrm{to}\quad & \boldsymbol{Ax} \leq \boldsymbol{b} \\
& \boldsymbol{x} \in \mathbb{R} \\
\end{align}
上記の問題例の場合、各パラメータは次のようになります。不等号の向きを変えるために、いくつかの制約条件には $-1$ をかけています。係数に注意してください。
線形計画問題(p.1)
\boldsymbol{c}^T = \begin{bmatrix} 100 & 300 \end{bmatrix},\quad
\boldsymbol{A} = \begin{bmatrix} -4 & -1 \\
-1 & -2 \\
-1 & -5 \\
0 & 1 \\
\end{bmatrix},\quad
\boldsymbol{b} = \begin{bmatrix} 22 \\ 13 \\ 29 \\ -1 \end{bmatrix}
これらのパラメータを scipy.optimize に与えれば、次のように求解することができます。もし、決定変数 $\boldsymbol{x}$ に非負条件があるときは x = [[0, None],[0, None]] のようにします。
scipy.optimizeによるLPの求解
import numpy as np
from scipy import optimize
c = np.array([100,300])
A = np.array([[-4,-1],[-1,-2],[-1,-5],[0,1]])
b = np.array([[22,13,29,-1]])
x = [[None, None],[None, None]]
opt = optimize.linprog(c, A, b, bounds=x)
print(f'x1* = {opt.x[0]:.2f}')
print(f'x2* = {opt.x[1]:.2f}')
print(f'f(x1*,x2*) = {opt.fun:.2f}')
実行結果
x1* = -2.33
x2* = -5.33
f(x1*,x2*) = -1833.33
決定変数の最適値 $x_1^{*}=-2.33$、 $x_2^{*}=-5.33$、目的関数の最小値 $f(x_1^{*},x_2^{*})=-1833.33$ が求まりました。
図法解法の解説のための作図
ここでは、横軸を $x_1$、縦軸を$x_2$ として作図するのですが、プログラムのなかでは、便宜上、横軸を $x$、縦軸を $y$ としています。実行可能領域の塗りには plt.fill_between を利用します。
図解法の作図(Python)
import numpy as np
import matplotlib.pyplot as plt
x1_ = -2.33 # 最適解（別途求めておく
x2_ = -5.33 # 最適解
xRange=(-10,2) # 作図範囲 X
yRange=(-10,2) # 作図範囲 Y
n = 200 # 作図分解能
x = np.linspace(xRange[0], xRange[1], n)
# 各式を x2 = a * x1 + b の形式で記述
of = -1/3 * x - 7.5 # 目的関数 切片は適当に
c1 = ( -4 * x - 22) / 1 # 制約条件 1
c2 = ( - x - 13) / 2 # 制約条件 2
c3 = ( - x - 29) / 5 # 制約条件 3
c4 = 0 * x - 1 # 制約条件 4
# 目的関数と制約条件を描画
plt.figure(figsize=(5,5), dpi=120)
plt.plot(x, c1 , color='0.0')
plt.plot(x, c2 , color='0.0')
plt.plot(x, c3 , color='0.0')
plt.plot(x, c4 , color='0.0')
plt.plot(x, of , color='red', linestyle='dashed')
plt.text(xRange[-1], of[-1], ' (0)', size = 8, verticalalignment='center')
plt.text(xRange[-1], c2[-1], ' (2)', size = 8, verticalalignment='center')
plt.text(xRange[-1], c3[-1], ' (3)', size = 8, verticalalignment='center')
plt.text(xRange[-1], c4[-1], ' (4)', size = 8, verticalalignment='center')
plt.hlines([0], xRange[0], xRange[1], linewidth=3)
plt.vlines([0], yRange[0], yRange[1], linewidth=3)
# 最適点をプロット
plt.plot(x1_, x2_, color='red', marker='o')
# 実行可能領域を塗る
yUb = c4
yLb = np.maximum(np.maximum(c1, c2), c3)
plt.fill_between(x, yUb, yLb, where=yUb>yLb, facecolor='b',alpha=0.3)
plt.grid(color='0.8') # グリッド
# 各軸のラベルとメモリ
plt.xlim(xRange[0], xRange[1])
plt.xticks(range(xRange[0], xRange[1]+1, 1))
plt.ylim(yRange[0], yRange[1])
plt.yticks(range(yRange[0], yRange[1]+1, 1))
plt.savefig('fig.svg',format = 'svg'); # SVG形式で出力
plt.show()
実行結果
ここに添付した画像は png ですが、実行フォルダ内に SVG 形式のファイルも出力されているので、あとは InkSpace などで微修正をします。