Pythonで線形計画問題を解く即戦力コード with Pyomo

はじめに

この記事では,

Python製のモデリングツールPyomoを使って標準的な線形計画問題を解く方法

を紹介します.
線形計画問題の数学的な話,解法アルゴリズム,ソルバー(gurobiとかcplexとか)の詳細な仕様については書きません.

Pythonで線形計画問題を解く方法は数多くあるのですが,その多さとPython界隈の開発のスピードのおかげで(せいで),参入障壁が高くなっているように感じます.

そもそもPythonで数値計算を行うことを目的としている人々にとっては,写経でもまず動くコードが欲しいのだと思います.
(自分もそうです.CSガチ勢でもなく,あくまでツールとして使いこなすことを目的にPythonを使っています,)

そこで,今回は自分が混乱したところを中心に,線形計画問題を解く即戦力コードを紹介したいと思います.

(*) 参考記事,サイトは数は最後にまとめます.

ソルバーとモデラー

最初の混乱ポイントががソルバーとモデラーの話です.

まず,ソルバーは汎用的な数理問題を解くツールです,
代表的なものに,Gurobi,cplex, ipoptなどがあります.

これらのソルバーは厳格な入力フォーマットを要求し,仕様も初心者に難しいので直接使うのは困難です.
さらに,大規模な問題でも高速に解くことを目的にしているため,PythonではなくCやC++で書かれています.

ここで,登場するのがモデラーです.モデリングツールとも呼ばれます.

これはインタプリンタ言語であるPythonとソルバーのつなぎ役を担います.

具体的には,ソルバーに渡すための入力フォーマットを作成し,ソルバーが解いてくれたらPythonの変数として受け取るなどの仕事をします.

これらは当然Pythonで書かれており,直感的なPython記法で記述が可能で,初心者にも一応優しい作りとなっています.

つまり,
モデラーを使ってPythonで問題(モデル)を記述し,モデラーがソルバーに投げて,解をモデラーが受け取る
というようなシステムのイメージです.

(*)文献によっては両者の違いを区別していない場合もあります,(明確な線引きは個別のツール依存で難しいため)

ソルバーの種類

正直,自分もあまり詳しくないのですが有名どころは以下です.
Scipy

  • 純Pythonのソルバー
  • 書きやすいが,大規模問題には無力(100変数ぐらいでストレスになる計算時間)

CBC

  • C#で書かれたソルバー
  • COIN-OR Projectの一貫
  • 混合整数計画問題に対する分枝切除法(branch and cut method)を実装している

glpk

ipopt

  • C++で主双対内点法を実装
  • 滑らか?な解がでる印象

Gurobi

  • 商用ソルバー,アカデミックライセンスあり
  • 選択可能なアルゴリズム Gurobi Method
    • 単体法
    • 内点法
    • バリアー法

cplex

モデラーの種類

線形計画問題を扱うPythonのモデラーも複数あります.代表的な3つと僕個人の印象をまとめます,
Scipy

  • Python製の最適化ツール(モデラー・ソルバー)
  • アルゴリズム部分(前述の区分で言うソルバー部分)もPythonなのでひたすら遅い
  • 単体法,内点法を選択可能
  • 書きやすいため,規模の小さい問題には適している

PuLP

  • Python製のモデラー
    • COIN-OR Projectの一貫
  • 呼出可能ソルバー
    • CBC
    • glpk
  • 2017年夏頃より開発が止まっているらしい
  • 日本語の解説本も存在するが,英語含めドキュメントは貧弱
  • myPuLPというPuLP自体のモデラーも存在しややこしい

Pyomo

  • Python製のモデラー
  • 呼出可能ソルバー
    • Gurobi
    • glpk
    • ipopt
    • cplex
  • Googleコミュニティが活発
  • ドキュメントは貧弱だがなんとか頑張れる

今回はPyomoを使った方法を紹介します

線形計画問題の基本形

以下に線形計画問題の基本形[P]を示します.
どんな線形計画問題もこの形に落とし込めるので,この解き方をマスターすればあとは微調整でなんとかなります.
[P]
$$
\max_{x} \qquad b^{\mathsf{T}} x
\\
\mbox{s.t.} \qquad A x \leq b
\\
\qquad \qquad x \geq 0
$$
ここで,

  • $x \in \mathbb{R}^{N}$ : 決定変数ベクトル
  • $b \in \mathbb{R}^{N}$ : 係数ベクトル
  • $A \in \mathbb{R}^{M \times N}$ : 係数行列
  • $c \in \mathbb{R}^{M}$ : 定数ベクトル

です.
今回は,より具体的に,$N=2, M=2$とし,
$$
b =
\begin{bmatrix}
2 \\
3
\end{bmatrix}
\\
A =
\begin{bmatrix}
4 & 3 \\
1 & 2
\end{bmatrix}
\\
c =
\begin{bmatrix}
120 \\
60
\end{bmatrix}
$$
として計算します.
(*) ちなみにこの問題は,和歌山大学のレジュメの2ページにあったチョコレート工場の生産問題です,

環境構築

Pyomoをinstallします.いずれかの方法でOKです.(*)Pyomo GitHub

  • pip install pyomo
  • conda install -c https://conda.anaconda.org/conda-forge pyomo

Pyomoは,各種ソルバーを自在に呼び出すことが可能ですが,ソルバーを使うためにはローカル環境にソルバーがinstallされている必要があります.

Gurobiやipoptは導入が少し手間なので,今回はgplkを使います.
glpkは以下のいずれかで導入可能です.

  • pip install glpk
  • conda install -c conda-forge glpk

実装

長くなりましたが,実装していきます.

In [1]:
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

import numpy as np
pyomoからglpkを呼び出せるか確認します.

In [2]:
# 使うソルバーの確認
SolverFactory("glpk").available() == True
Out[2]:
True
無事,使えることがわかったので,さきほどの線形計画問題[P]を作っていきます.

In [3]:
#変数・パラメータの設定

#係数行列の次元
N=2
M=2

#係数ベクトル,係数行列の作成
b = np.array([2,3], dtype="float")
A = np.array([[4,3],[1,2]], dtype="float")
c = np.array([120,60], dtype="float")
ここからが曲者です.Pyomoでは決定変数にもPythonの変数を用意します.
そして謎の初期化関数もセットで用意します,(理由はよくわかりません)

In [4]:
#決定変数の設定
x = np.zeros(shape=(N), dtype="float")

#決定変数初期化用関数---------------------
xx = {i:x[i] for i in range(len(x))}
def init_x(model,i):
    return xx[i]
# -----------------------------------
さらに,謎のindex_listも用意します.本当に単純なリストです.

In [5]:
#index list
N_index_list = list(range(N))
M_index_list = list(range(M))
不等式制約と,目的関数を以下のように記述します.
ここら辺はこう書くものと割り切っています.ここで.index_listが活躍します.
(おそらくPyomoは一般的な配列のindexを扱えないのだと思います,そのため単純なindex_listによって要素指定をします.)

In [6]:
#不等式制約
def const(model, j):
    tmp = sum(A[j,i]*model.x[i] for i in N_index_list)
    return tmp <= c[j]

#目的関数
def func(model):
    return sum(b[i] * model.x[i] for i in N_index_list)   
Pyomoのモデルに制約や決定変数を組み込む流れは以下の通りです.

In [7]:
#モデル
model = pyo.AbstractModel("LP-sample")

#決定変数をモデルに組み込み
model.x = pyo.Var(N_index_list, domain=pyo.NonNegativeReals, initialize=init_x)

#目的関数をモデルに組み込む
model.value = pyo.Objective(rule=func, sense=pyo.maximize)

#制約をモデルに組み込み
model.const = pyo.Constraint(M_index_list, rule=const)
ソルバーをgplkに設定します,ソルバーオプション(適用アルゴリズムなど)についても詳細に設定できますが,その方法は今回は割愛します.
デフォルトでもいい感じに解いてくれます.

In [8]:
#ソルバーの選択とオプション
solver_name = "glpk"
opt = SolverFactory(solver_name)
さらに以下の手順で問題をソルバーに投げて,解を受け取ります.

In [9]:
instance =  model.create_instance()
instance.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) #双対変数取得用に必要
results = opt.solve(instance, tee=True)#tee=Trueでログ出力
instance.load(results)
results.write()
GLPSOL: GLPK LP/MIP Solver, v4.65
GLPK Simplex Optimizer, v4.65
3 rows, 3 columns, 5 non-zeros
Preprocessing...
2 rows, 2 columns, 4 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  4.000e+00  ratio =  4.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 2
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (2)
*     2: obj =   9.600000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (40412 bytes)


# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 96.0
  Upper bound: 96.0
  Number of objectives: 1
  Number of constraints: 3
  Number of variables: 3
  Number of nonzeros: 5
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.03223729133605957
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
本物の出力はもっと長いかもしれません.(本記事ではディレクトリの場所がわかってしまう出力を削除しています)

In [10]:
#解の概要を見る
instance.display()
Model LP-sample

  Variables:
    x : Size=2, Index=x_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  12.0 :  None : False : False : NonNegativeReals
          1 :     0 :  24.0 :  None : False : False : NonNegativeReals

  Objectives:
    value : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True :  96.0

  Constraints:
    const : Size=2
        Key : Lower : Body  : Upper
          0 :  None : 120.0 : 120.0
          1 :  None :  60.0 :  60.0
In [11]:
#最適値をだけ見る
print(instance.x[0].value)
print(instance.x[1].value)
12.0
24.0
In [12]:
#最適目的関数値を見る
instance.value.expr()
Out[12]:
96.0
In [13]:
#双対変数の一覧を見る
instance.dual.display()
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
    Key      : Value
    const[0] :   0.2
    const[1] :   1.2
In [14]:
#双対変数値だけ見る
print(instance.dual[instance.const[1]])
print(instance.dual[instance.const[0]])
1.2
0.2

おわりに

おおむね即戦力になるコードの紹介ができたと思います.

ただ,Pyomoもドキュメントが少し怪しかったり,美しい唯一の書き方がなかったりするので自分なりに試行錯誤は必要です.

違うソルバーの導入方法や使い方,オプションの設定方法についてはまた記事にしたいと思います,

参考サイト

0