しばらく，線形計画でできる話をしていきます。

線形計画問題とは

線形「計画」とはいいますが Linear Programming の訳語です。リニアが線形（一次式）で，プログラミングはこういう文脈では「計算」くらいの意味合いです。どんな計算かといえば，

何本かの一次不等式で表される領域のなかで
一次式の目的関数の値を最大化 (or 最小化) するような変数の組を特定する

計算であり，要するに最大最小問題の一種です。ツウの人たちはLPっていいますね。

大学入試でも，2変数の線形計画は出題されることがありますね！

mathtrain.jp

大学入試なんかだと，「最大値を求めよ」と出題され，そのときの変数の値はオマケ的だけど必ず答えにつけようねとお作法みたいに思われていますが，実用的には，「最適のときの変数の値がコレです！！」というのをちゃんと求めるのが超重要です。

ちゃんと書きますと， $n$ 個の変数 $\displaystyle{ x_i (i = 1, \cdots, n)}$ に対して， $m$ 本の制約条件（一次不等式）で表される領域 $D:$

$\displaystyle{\sum_{i = 1}^{n} a _ {ij} x _ i \le b_j \quad(j=1, \cdots, m)}$

の中で，目的関数

$\displaystyle{z = \sum_{i = 1}^{n} c_i x_i}$

を最大（または最小）とするのが線形計画問題です。

以上をまとめて，

$\begin{aligned} \text{maximize} \quad &z = \sum_{i = 1}^{n} c_i x_i \\ \text{subject to} \quad &\sum _ {i = 1}^{n} a _ {ij} x _ i \le b _ j \quad (j=1, \cdots, m) \end{aligned}$

のように書きます。

最大値と最小値の両方を求めなくていいの？と思う方もいるでしょう。が，実際に問題を組んでみると，最大値が $+\infty$ だったり，ただの0だったりすることも多いので，意味のあるほうだけ求めるのが通常です。もちろん，必要があれば両方求めてOKです。

ちなみに，制約条件に等式制約を使うのもOKです。（スラック変数というものを持ち出すと不等式制約をすべて等式制約に変換できます。逆に，等式制約 $A=B$ も不等式制約 $A \le B \land A \ge B$ に変換できますが細かい話は割愛。どっちも使えます。）（鉄道方面から来た方，曲線に付けるカントとスラックのスラックと一緒です）

線形計画問題の例題を解く

今回は2変数の簡単な例題を，PuLPで解いてみましょう。

PuLPについて

PuLPは，Pythonで線形計画（と整数計画）に使えるモジュールです。セットアップについては他に詳しいサイトがたくさんありますので今回は割愛。

個人的には，Anacondaを使わずに，VS Code内で完結するようなセットアップが好みです。とりあえずはpipでjupyterとpulpを入れてやればいいと思います。

例題

$\begin{aligned} \text{maximize} \quad &z = 4x + 5y \\ \text{subject to} \quad & x + y \le 3 \\ & x + 2y \le 4 \\ & x \ge 0, \quad y \ge 0 \end{aligned}$

さきほどの「高校数学の美しい物語」さんのページに載っている例題と一緒です。2変数なので図形的にも解け，最適解は $(2,1)$ で $z=13$ です。

PuLPの使い方

PuLPを使うときのおおまかな流れは

インポートする (from pulp import *)
モデルオブジェクトを宣言する
変数オブジェクトを宣言する
目的関数をセットする
制約条件をセットする
最適化する
結果を見る

です。

手順1,2,3は，はじめはおまじないだと思っても結構です。

from pulp import *

# モデルオブジェクトの宣言
model = LpProblem("sample", sense=LpMaximize)

# 変数オブジェクトの宣言
x = LpVariable("x", lowBound=0)
y = LpVariable("y", lowBound=0)

モデルオブジェクトで，sense=LpMaximizeは最大化問題であることを指定しています。デフォルトだと最小化問題になりますが，最大化問題のつもりで最小化問題を解いたりするポカがあるのでいつも書くようにしておくといいです。

変数オブジェクトで，lowBound=0は変数の下界を0にします。（デフォルト $-\infty$ です。）変数の上界がある場合はupBound=で設定できますがやはりデフォルトは $+\infty$ です。

手順4と5で目的関数と制約条件をセットします。同じ書き方をしますが順番が大事です。

# 最初に+=したのが目的関数と解釈される
model += 4*x + 5*y, 'objective'

# 2番目以降は制約条件と解釈される
model += x +   y <= 3, 'c1'
model += x + 2*y <= 4, 'c2'

objective, c1, c2 は式につけた名前です。

ここで一度modelの中身をprintしてみたりなんかすると，中身が見えます。

sample:
MAXIMIZE
4*x + 5*y + 0
SUBJECT TO
c1: x + y <= 3

c2: x + 2 y <= 4

VARIABLES
x Continuous
y Continuous

この問題で合っていれば，手順6 最適化してあげます。

model.solve()
print("Status:", model.status, LpStatus[model.status])

.solve()したあと，status属性を確認してあげて，最適に解けたかどうか確かめます。model.status == 1なら最適解が見つかったという意味です。LpStatusは辞書なのでキーにstatusを入れてあげると意味の取れる文字列で見せてくれます。

Status: 1 Optimal

が表示されたのでOKでした。最適値があるはずなのに Infeasible (実行不可能) と言われたり Unbounded (無限解) と言われたら，たぶん定式化（目的関数と制約条件の設定のこと）でミスっています。

手順7 最大値とそのときの各変数の値を確認してあげます。

print("z =", value(model.objective))
for var in model.variables():
    print(var, var.varValue)

value()関数を噛ませたり，varValue属性を読んだりするのがややこしく感じるかもしれません。変数オブジェクトはいつものfloatオブジェクトとはモノが異なるのでナマで参照しても値のようには振る舞いません。

z = 13.0
x 2.0
y 1.0

と，図形的に解いたのと同じ答えが求められました。

コードの全文をまとめると以下の通り。

from pulp import *

# モデルオブジェクトの宣言
model = LpProblem("sample", sense=LpMaximize)

# 変数オブジェクトの宣言
x = LpVariable("x", lowBound=0)
y = LpVariable("y", lowBound=0)

# 最初に+=したのが目的関数と解釈される
model += 4*x + 5*y, 'objective'

# 2番目以降は制約条件と解釈される
model += x +   y <= 3, 'c1'
model += x + 2*y <= 4, 'c2'

# 求解とその前後処理
print(model)
model.solve()
print("Status:", model.status, LpStatus[model.status])

# 結果の確認
print("z =", value(model.objective))
for var in model.variables():
    print(var, var.varValue)

線形計画とは結局なんなのか

線形計画問題 (LP) がもてはやされるのはなぜか。

$n$ 変数のLPでも高速に解けるアルゴリズムが存在する。
工学的な応用範囲が広い。

もちろん，「残念ながら，世の中は非線形でできている」のはみんな知っているのですが，線形で表せる工学的（工業的）な営みや，近似して線形だと思えればそれでいいだろうというテーマがたくさんあります。組合せ的な問題も，なかにはLPになったりするものもあるのです。

LPと同じ形で、変数が整数値しか取れない場合のものを整数計画問題といい、これもPuLP (に同梱のCBCソルバー) で解くことができます。ナップサック問題のような組合せ最適化問題が整数計画に帰着されます。k-means、k-centerのようなクラスタリングや、みんな大好きTSPも整数計画になります。（ところが整数計画シリーズは、「整数計画問題にできても、現実的な時間で最適解を求められない」ことがあります。）

Python + PuLP での最適化は、並木誠先生の本があり、理論の説明から入っていけますがミスプリが多いというAmazonのレビューもあります…。

Pythonによる数理最適化入門 (実践Pythonライブラリー)

作者:並木誠
出版社/メーカー: 朝倉書店
発売日: 2018/04/09
メディア: 単行本

正誤表とサンプルコードは出版社のHPでゲットできます。

補遺

フリーで使えるソルバーということでPuLP（に同梱のCBCソルバー）を当分使っていきますが，私は実はGurobi使いです（いまのところ）。pythonAPIの出来が良く、問題の記述がPuLPより簡単で柔軟で分かりやすいだけでなく，求解の（速度的な）性能も滅茶苦茶に強いです。アカデミアにいる人はアカデミックライセンスがもらえるので公式サイトを訪問してください。

python API の使い勝手でいうとSCIPもGurobiとかなり似てて使いやすそうです。が，ライセンスがハッキリしません（アカデミアは無料，商用は連絡してね（有料？），なんだけど個人の趣味用途は？？）。

CBCはその名を Coin Branch and Cut Solver というだけあって，本当は分枝カット法に基づく高度な解法を実装できるんじゃないかと思うんですが，残念ながらPuLPはそういうところまで踏み込めないみたいです。C++が使える人はゴリゴリにCBCを使い倒してみてもいいんじゃないでしょうか。pythonで比較的深いところまでCBCを叩けるcbcpyなるモジュールもあるらしいので，このブログでもそのうちやってみたいことの一つです。

Googleのor-toolsも無料にしては優秀なんじゃないかと思っていますが，Python APIがC++のコードから自動生成したもののようでエラく書きにくいので，今回はパス。が，便利そうなツールなのでこれもそのうちやってみたいですね。（ただの混合整数計画を組んでしまうと結局回るのはCBCソルバーだそうです。）こちらもC++が読み書きできる方にはよいかもしれません。

Local Optima

数理最適化を趣味的に楽しむ。

Python + PuLPで線形計画

線形計画問題とは

線形計画問題の例題を解く

PuLPについて

例題

PuLPの使い方

線形計画とは結局なんなのか

補遺