線形計画法の勉強メモ

私は, 現在工学系大学院で最適化の講義を履修しています. 
講義の自習として勉強メモを記述します. 

本記事は, 線形計画法の勉強メモです. 

## 線形計画問題とは
以下は, 「しっかり学ぶ数理最適化」[1]からの引用となります. 
> 線形計画問題は目的関数が線形関数で, すべての制約条件が線形の等式もしくは不等式で表された最適化問題である. 

また, 「Pythonではじめる数理最適化」[2]では, 次のようにも紹介されています. 
> 線形計画問題：変数が実数値をとる問題。生産量（0.3kgなどを許す）を決める場合など。

## 例題
ここでは, 次のような例題を考えます. 
> $x, y$が4つの不等式 $x\geq0, \quad y\geq0, \quad x-3y+5\geq0, \quad 2x+5y-27\leq0$を満たすとき, $x+2y$のとる値の最大値を求めよ. 

線形計画問題の定義通り, 目的関数$x+2y$は線形関数となっています. 
また, 制約条件 $x\geq0, \quad y\geq0, \quad x-3y+5\geq0, \quad 2x+5y-27\leq0$ は線形の不等式で表されていることがわかります. 

この例題について, 高校数学まででの解法とPythonのPuLPライブラリそれぞれを使って解を求めたいと思います. 

## 解 ( 高校数学まででの場合 )
4つの不等式 $x\geq0, \quad y\geq0, \quad x-3y+5\geq0, \quad 2x+5y-27\leq0$によって示される領域$D$は, 次の図の通りとなる. 

![Plot1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3582223/7c9abb53-3d6b-2658-16c4-b017c6ec110a.png)
( 画像はWolframeAlphaでプロットしました. )

ここで, $x+2y=k$とすると, これは傾き$-\frac{1}{2}$, $y$切片$\frac{k}{2}$の直線を表す. 
この直線と領域$D$が共有点を持つような$k$の値において, 最大となるようなものを求めることを考える. 
これは, 直線の傾き$-\frac{1}{2}$より, 点$\biggl( \frac{27}{2}, 0 \biggr)$を通ると分かる. 
よって, $x+2y$は, $x=\frac{27}{2}, y=0$のとき, 最大値$\frac{27}{2}$をとる. 

![Plot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3582223/7f1951bb-a3aa-4f48-844a-90645e0bf483.png)
( これも画像はWolframeAlphaでプロットしました. 画像がちょっと見にくい... )

## 解 ( PythonのPuLPライブラリを用いた場合 )
続いて, PythonのPuLPライブラリを用いて解を求めていきます. 

今回は, Google Colaboratory上で実行しました. 

まずは, Pythonのバージョンを確認していきます. 
```python
$ !python --version
> Python 3.10.12
```

続いて, PuLPライブラリをインストールします. PuLPライブラリはColaboratoryにデフォルトで入っているわけではないみたいで, インストールせずに実行するとエラーが発生します. 
```python
$ !pip install pulp
```
適切にインストールされていれば大丈夫です. 多分問題ないはず...

それでは, いよいよ線形計画問題を解いていきます. コードは「Pythonではじめる数理最適化」[2]を参考にしました. 
```python
import pulp

prob = pulp.LpProblem('LinearProgramming', pulp.LpMaximize)

x = pulp.LpVariable('x', cat='Continuous')
y = pulp.LpVariable('y', cat='Continuous')

prob += x - 3*y + 5 >= 0
prob += 2*x + 5*y - 27 <= 0
prob += x >= 0
prob += y >= 0
prob += x + 2*y

status = prob.solve()

print(pulp.LpStatus[status])
print('x: ', x.value(), ', y: ', y.value(), ', max: ', prob.objective.value())
```

1行目でPuLPライブラリをインポートしています.

3行目では, 問題を設定するための数理モデルを定義しています. 今回は最大化問題のため, LpProblemの第2引数にpulp.LpMaximizeを設定しています (最小化問題の場合は, pulp.LpMinimize). また, 第1引数に指定する名前は任意のもので良いみたいなので, 今回は線形計画問題を示すLinearProgrammingとしています. 

5,6行目では, 今回用いる変数$x$と$y$の定義を行っています. ここも, pulp.LpVariableの第1引数は任意の名前みたいです. また, 今回は連続な関数を用いることを考えているため, 第2引数にContinuousを設定しています. 

8~11行目では, 制約条件を定義しています. 今回で言えば, 領域$D$を形成する条件達ですね. problem += ( 右辺 )の形になっていて, この( 右辺 )の部分に素直に制約条件を書いていけます.

12行目では, 目的関数を定義しています. 制約条件との違いは, ( 右辺 )部分に$=$や$\geq$などがないことですね. 

14行目では, いよいよ線形計画問題を解いていきます. 

16, 17行目では, 結果を記述しています. pulp.LpStatus[status]は数理モデルが解けたかどうかを示しており, 17行目で実際の解 (ここでは, 最大値とそのときの座標) を示しています. 


このコードを実行した結果は, 次の通りです. 
```python
Optimal
x:  13.5 , y:  0.0 , max:  13.5
```

## 所感
ということで, どちらの解も同じとなりました. 2変数であれば高校数学までの解法でもできますが, それ以上だとやはりPythonで実行した方が楽そうですね. 

## 参考文献
[1] 梅谷俊治, "[しっかり学ぶ数理最適化 モデルからアルゴリズムまで](https://www.kspub.co.jp/book/detail/5212707.html)", 講談社, 2020.10.23, p.13.
[2] 岩永次郎, 石原響太, 西村直樹, 田中一樹, "[Pythonではじめる数理最適化 -ケーススタディでモデリングのスキルを身につけよう-](https://www.ohmsha.co.jp/book/9784274227356/)", オーム社, 2021.9.20, pp.16-21.