線形計画問題の図解法を解説するための資料をmatplotlibで作成（Python）

概要

　線形計画問題の図解法を解説するための資料を 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 などで微修正をします。