はじめに

本記事の内容

Pythonで数理計画問題（最適化問題）を解く方法については、既に様々な記事が出ています。ただその多くは数個の変数と数個の制約式からなるシンプルな例題を対象としており、変数や制約式が添え字で展開されるような"実用的"な規模の問題に適用するには、そのままではわかりにくいことがあります。

本記事では、こういった"実用的"な変数・制約式から成る数理計画問題を対象として、Pythonで記述を行う方法を示します。このとき、複数のモデル記述言語（モデラー）を比較することで、それぞれの特徴を示します。

本記事に含まれないこと

「Pythonで数理計画問題を解く」といっても、例えば最急降下法や遺伝アルゴリズム等の解法実装の実装方法を示しているわけではありません。「モデラー」を用いて数理計画問題の目的関数と制約式を記述し、その解を求めるのは汎用の「ソルバー」に任せるという方法をとっています。
また、数理計画問題についての解説は特に行っていません。数理計画問題を、通常の数学記法で行えていることは前提とします。

環境

以下は、Google Colaboratory上でのコーディングを示しています。

モデル記述言語（ライブラリ）

モデル記述のためのライブラリとして以下を取り上げます。これら以外にもCVXOPT、PICOSといったライブラリもあります。
世の中にはPuLPの利用例が多いようですが、モデルに整数変数や非線形項を入れる可能性があるのであれば、最初からpyomoを選択することがおすすめです。

ライブラリ	解ける問題形式	参考資料
PuLP	LP	記述法¹ シャドウプライスの取り方²
PYTHON-MIP	LP, MIP	記述法³ シャドウプライスの取り方⁴
pyomo	LP, MIP, MINLP, NLP ,QCP, MIQCP	ソルバーインストール方法⁵ シャドウプライスの取り方⁶

例題

本記事では、以下の線型の数理計画問題を例として、各モデル記述言語（モデラー）で記述していきます。
なお、変数は大文字、係数や添え字は小文字で記載しています。

電力供給問題

各時刻$t$において、複数の発電所$g$を用いて、総発電コスト$COST$が最小になるよう電力を供給したい。各発電所からどのように発電量$GEN$を出力すればよいか。
ただし、各時刻の発電量合計は需要$demand$に一致する必要がある。また各発電所の発電量には上限$cap$があり、また発電量の時刻間変化率には係数$coeff$の上下限がある。

$ COST \rightarrow min $

$ subject \ to $
$\quad COST = \sum_{g=0}^{num_{gens}} \sum_{t=0}^{num_{hours}} cost_g*GEN_{g,t} $
$\quad \sum_{g=0}^{num_{gens}} GEN_{g,t} = demand_t \quad (\forall\ t)$
$\quad GEN_{g,t} \leqq cap_g \quad (\forall\ g,t)$
$\quad GEN_{g,t} \leqq GEN_{g, t-1} \times (1+coeff_g,t) \quad (\forall\ g,\quad t\geqq1)$
$\quad GEN_{g,t} \geqq GEN_{g, t-1} \div (1+coeff_g,t) \quad (\forall\ g,\quad t\geqq1)$

記述の比較

共通部分

まず、モデラーによらない共通部分。

共通ライブラリ

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

数値例

ここでは、3つの発電所（発電容量はいずれも10）で、需要（最小10、最大20）を満たすように電力を供給するという数値例を用います。（最もコストの高い発電所は不要だと思うでしょうか？）
これらのパラメータはスカラーあるいはリストの形式で持っておきます。

num_gens = 3
num_hours = 6
demand = [12, 10, 18, 20, 14, 12]
capacity = [10, 10, 10]
cost = [5, 10, 15]
coeff = [0, 0.5, 2]

イテレータ

モデル記述に使用するイテレータを定義しておきます。

generators = range(num_gens)
hours = range(num_hours)
hours1 = range(1, num_hours)

ライブラリの読み込み

以降は、各モデラーでの記述方法を、並列して記載しています。

PuLP

!pip install pulp
from pulp import *

PYTHON-MIP

!pip install mip
from mip import *

pyomo

ソルバー（CBC）もここでインストールしています。

!pip install pyomo
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
!apt-get install -y -qq coinor-cbc

モデルの生成

いずれも、モデルの空箱のようなものを作り、そこに変数や制約式を追加していくという方法をとります。
まずはモデルの空箱作りです。

PuLP

m1 = pulp.LpProblem()

PYTHON-MIP

m1 = mip.Model()

pyomo

m1 = pyo.ConcreteModel()

なおpyomoには、「具体モデル」(ConcreteModel)と、「抽象モデル」(AbstractModel)という概念があります。抽象モデルでは、モデル記述中には式の具体的な係数の内容（添え字の上限や値そのもの）を与えず、解くときに別途外部のデータファイルを読み込みます。ただし具体モデルでもある程度似たことはできるので、ここでは具体モデルでの記述を行っています。

変数の生成

ここからそれぞれの違いが際立ってきます。特に、添え字付き変数$GEN_{g,t}$の表現に着目してください。

PuLP

PuLPは、いきなりここが最も複雑です。
添え字を文字列として組み込んだ変数「GEN_0_0」「GEN_0_1」…を個別に生成します。
このとき、まず添え字の全組み合わせを生成するため、Python関数itertools.product()を使用しています。これをイテレータとしてforで繰り返しながら、演算子%（C言語のsprintf関数のような動作）を用いて文字列の指定箇所にg,tの値を代入した変数名を持つ変数を作り、かつ添え字(g, t)とその変数の対応を示した辞書(dict型)変数を作っています。

GEN = {(g,t):pulp.LpVariable('GEN_%d_%d'%(g,t)) for g,t in list(itertools.product(generators, hours))}
COST = pulp.LpVariable('COST')

なおGENの中身を見てみると、次のような辞書形式になっていることがわかります（変数が文字列で表現されていますが、実際はLpVariableです）。

{(0, 0): GEN_0_0, (0, 1): GEN_0_1, (0, 2): GEN_0_2, (0, 3): GEN_0_3, (0, 4): GEN_0_4, ...

PYTHON-MIP

add_var_tensorを使うことで、添え字付き変数をシンプルに生成することができます。

GEN = m1.add_var_tensor((num_gens, num_hours), 'GEN')
COST = m1.add_var('COST')

pyomo

添え字付き変数とそうでない変数を、同一の関数Varでシンプルに生成することができます。

m1.GEN = pyo.Var(generators, hours)
m1.COST = pyo.Var()

目的関数の追加

PuLP

最適化の方向を定めたあと、+=の記法でモデルに目的関数や制約式を追加していきます。等式・不等式ではない（すなわち等号・不等号を含まない）数式を追加すると、自動的に目的関数であると認識されます。
+=記法がなんとなく気持ち悪ければ、objectiveプロパティへの代入も可能です。

m1.sense = LpMinimize
m1 += COST
# 以下も同じ意味
# m1.objective = COST

PYTHON-MIP

最適化の方向を定め、目的関数を定義します。
なお、後述のとおり制約式はPuLPと同様+=での追加ができるのですが、目的関数は明示的にobjectiveプロパティへ代入する必要があります。

m1.sense = mip.MINIMIZE
m1.objective = COST
# 以下も同じ意味
# m1.objective = mip.minimize(COST)

pyomo

pyomoの記述方法は特徴的です。
まず数式を返す関数を定義し、次にその関数をObjective関数により目的関数としてモデルへ登録する、という2段階の手順を踏みます。

def obj(m):
  return m.COST
m1.obj = pyo.Objective(rule=obj, sense=pyo.minimize)

制約式の追加１

まず、COSTを計算するための数式の記述方法を比較します。
$ COST = \sum_{g=0}^{num_{gens}} \sum_{t=0}^{num_{hours}} cost_g*GEN_{g,t} $

PuLP

Pythonに特徴的な「リスト内包表記」が使われていますが、数式どおりの記述ともいえます。
カンマ以下は制約式の名前で、省略することもできますが、最適解を求めた後に数式を参照する際（シャドウプライスやスラック変数の取得等）には名前があるほうが便利です。
+=の記法ではなく、addConstraint関数で追加していくこともできます。

m1 += COST == pulp.lpSum(pulp.lpSum(cost[g] * GEN[g,t] for g in generators) for t in hours), 'totalcost'
# 以下でも同じ意味
# m1.addConstraint(COST == pulp.lpSum(pulp.lpSum(cost[g] * GEN[g,t] for g in generators) for t in hours), name='totalcost')

PYTHON-MIP

関数名がlpSumでなくxsumである点を除けば、PuLPと同じです。
+=の記法ではなく、add_constr関数で追加していくこともできます。

m1 += COST == mip.xsum(mip.xsum(cost[g] * GEN[g,t] for g in generators) for t in hours), 'totalcost'
# 以下でも同じ意味
# m1.add_constr(COST == mip.xsum(mip.xsum(cost[g] * GEN[g,t] for g in generators) for t in hours), name='totalcost')

pyomo

目的関数と同じく、pyomoの特徴的な記述となっています。等式や不等式を返す関数を定義し、その関数をConstraint関数により制約式としてモデルへ登録します。
関数の等式の中の記述自体は、PuLPやPYTHON-MIPと似ています（なおこのsumはpyomoライブラリの関数ではなく、一般のsumです）。
制約式の名前は必須です（m1.totalcostの「totalcost」が制約式の名前に当たります。今回は関数名totalcostと同じにしています）。

def totalcost(m):
  return m.COST == sum(sum(cost[g] * m.GEN[g,t] for g in generators) for t in hours)
m1.totalcost = pyo.Constraint(rule=totalcost)

制約式の追加２

残りの数式4種を記述します。COSTの計算式と異なり、それぞれ複数の添え字に対して追加すべき数式です。

$ \sum_{g=0}^{num_{gens}} GEN_{g,t} = demand_t \quad (\forall\ t)$
$ GEN_{g,t} \leqq cap_g \quad (\forall\ g,t)$
$ GEN_{g,t} \leqq GEN_{g, t-1} \times (1+coeff_g,t) \quad (\forall\ g,\quad t\geqq1)$
$ GEN_{g,t} \geqq GEN_{g, t-1} \div (1+coeff_g,t) \quad (\forall\ g,\quad t\geqq1)$

PuLP

添え字の組み合わせを、明示的にfor文で繰り返します。
制約式の名前も繰り返されるたびに生成します。ここでも演算子%を使っています。

for t in hours:
    m1 += pulp.lpSum(GEN[g, t] for g in generators) == demand[t], 'demand_balance_%d'%t

for g in generators:
  for t in hours:
    m1 += GEN[g, t] <= capacity[g], 'capacity_constraint_%d_%d'%(g,t)

for g in generators:
  for t in hours1:
    m1 += GEN[g, t] <= GEN[g, t-1]*(1+coeff[g]), 'load_following_1_%d_%d'%(g,t)
    m1 += GEN[g, t] >= GEN[g, t-1]*(1/(1+coeff[g])), 'load_following_2_%d_%d'%(g,t) #除算記号だけではエラー、分数を乗じる

PYTHON-MIP

PuLPとほぼ同じです。なお、明示的にfor文を書かなくてもコロンで代用できる（pandas等と同様の記述）場合があるので、一つ目の制約式は略して書いてみています。

for t in hours: #略して書ける
    m1 += mip.xsum(GEN[:, t]) == demand[t], 'demand_balance_%d'%t

for g in generators:
  for t in range(1, num_hours):
    m1 += GEN[g, t] <= capacity[g], 'capacity_constraint_%d_%d'%(g,t)

for g in generators:
  for t in hours1:
    m1 += GEN[g, t] <= GEN[g, t-1]*(1+coeff[g]), 'load_following_1_%d_%d'%(g,t)
    m1 += GEN[g, t] >= GEN[g, t-1]/(1+coeff[g]), 'load_following_2_%d_%d'%(g,t)

pyomo

等式や不等式を返す関数を定義し、その関数を制約式としてモデルへ登録します。
定義する関数の第2引数以降では添え字の各要素を受け取れるようにしておき、Constraint関数で登録する際に、添え字の範囲を表すイテレータを順番に指定します。for文として繰り返しの記述を書くわけではありません。
また、添え字付きの制約式名を明示的に生成する必要はありません。

def demand_balance(m, t):
  return sum(m.GEN[g, t] for g in generators) == demand[t]
m1.demand_balance = pyo.Constraint(hours, rule=demand_balance)

def capacity_constraint(m, g, t):
  return m.GEN[g, t] <= capacity[g]
m1.capacity_constraint = pyo.Constraint(generators, hours, rule=capacity_constraint)

def load_following_1(m, g, t):
  return m.GEN[g, t] <= m.GEN[g, t-1]*(1+coeff[g])
m1.load_following_1 = pyo.Constraint(generators, hours1, rule=load_following_1)

def load_following_2(m, g, t):
  return m.GEN[g, t] >= m.GEN[g, t-1]/(1+coeff[g])
m1.load_following_2 = pyo.Constraint(generators, hours1, rule=load_following_2)

この2段階のpyomoの記述方法は無駄なようにも思えますが、異なる制約式群から成る複数のモデルを並列で作成するときに、関数が再利用できるのです。

#（参考）制約式の1つ少ない別モデルを作ってみる
m2 = pyo.ConcreteModel()

#変数の宣言
m2.GEN = pyo.Var(generators, hours)
m2.COST = pyo.Var()

#目的関数と制約式の設定　ruleの関数は宣言済み
m2.obj = pyo.Objective(rule=obj, sense=pyo.minimize)
m2.totalcost = pyo.Constraint(rule=totalcost)
m2.demand_balance = pyo.Constraint(hours, rule=demand_balance)
m2.capacity_constraint = pyo.Constraint(generators, hours, rule=capacity_constraint)
m2.load_following_1 = pyo.Constraint(generators, hours1, rule=load_following_1)

ソルバーの指定と求解

各モデラーが対応しているソルバーは少しずつ異なります。
ここでは共通でCBCを使っていますが、この程度の規模のモデルであれば、何を使っても瞬く間に同じ解に収束するでしょう。

PuLP

solver = pulp.PULP_CBC_CMD()
res = m1.solve(solver)

PYTHON-MIP

m1.solver_name='CBC'
res = m1.optimize()

pyomo

後でシャドウプライスを取得するためには事前準備が必要のようです。

#後でシャドウプライスにアクセスできるようにしておく
m1.dual = pyo.Suffix(direction = pyo.Suffix.IMPORT)

solver = pyo.SolverFactory('cbc', executable='/usr/bin/cbc')
res = solver.solve(m1)

解へのアクセス１

多数の変数・制約式から成る問題は、ソルバーの返り値をそのままprintしてみて必要な解を目で探すのは非効率です。
まず、全変数とその値、全制約式とそのシャドウプライス（限界費用、潜在価格ともいう）にアクセスするための方法を示します。

PuLP

print('変数とその値')
for v in m1.variables():
  print(v.name, '=', v.varValue)

print('制約式とそのシャドウプライス')
for name, c in m1.constraints.items():
  print(name, ':', c, ':', c.pi)

出力例

変数とその値
COST = 898.88889
GEN_0_0 = 4.0
GEN_0_1 = 4.0
（略）
制約式とそのシャドウプライス
totalcost : COST - 5*GEN_0_0 - 5*GEN_0_1 - 5*GEN_0_2 - 5*GEN_0_3 - 5*GEN_0_4 - 5*GEN_0_5 - 10*GEN_1_0 - 10*GEN_1_1 - 10*GEN_1_2 - 10*GEN_1_3 - 10*GEN_1_4 - 10*GEN_1_5 - 15*GEN_2_0 - 15*GEN_2_1 - 15*GEN_2_2 - 15*GEN_2_3 - 15*GEN_2_4 - 15*GEN_2_5 = 0 : 1.0
demand_balance_0 : GEN_0_0 + GEN_1_0 + GEN_2_0 = 12 : 15.0
（略）

なお、制約式に名前を付けなかった場合、式には「_C2」のように追加したときの連番付きの名称が自動でつきます。

PYTHON-MIP

print('変数とその値')
for v in m1.vars:
  print(v.name, '=', v.x)

print('制約式とそのシャドウプライス')
for c in m1.constrs:
  print(c, ':', c.pi)

出力例

変数とその値
GEN_0_0 = 4.0
GEN_0_1 = 4.000000000000001
（略）
COST = 898.8888888888888
制約式とそのシャドウプライス
totalcost: -5.0 GEN_0_0 -5.0 GEN_0_1 -5.0 GEN_0_2 -5.0 GEN_0_3 -5.0 GEN_0_4 -5.0 GEN_0_5
	 -10.0 GEN_1_0 -10.0 GEN_1_1 -10.0 GEN_1_2 -10.0 GEN_1_3 -10.0 GEN_1_4 -10.0 GEN_1_5
	 -15.0 GEN_2_0 -15.0 GEN_2_1 -15.0 GEN_2_2 -15.0 GEN_2_3 -15.0 GEN_2_4 -15.0 GEN_2_5
	 +1.0 COST = -0.0 : 1.0
demand_balance_0: +1.0 GEN_0_0 +1.0 GEN_1_0 +1.0 GEN_2_0 = 12.0 : 15.0
（略）

なお、制約式に名前を付けなかった場合、式には「constr(1)」のように追加したときの連番付きの名称が自動でつきます（PuLPの制約式連番と異なり0開始）。

pyomo

print('変数とその値')
for v in m1.component_data_objects(ctype=pyo.Var):
    print(v.name, '=', v.value)

print('制約式とそのシャドウプライス')
for c in m1.component_data_objects(ctype=pyo.Constraint):
    print(c.name, ':', c.body, ':', m1.dual[c])

出力例

変数とその値
GEN[0,0] = 4.0
GEN[0,1] = 4.0
（略）
COST = 898.88889
制約式とそのシャドウプライス
totalcost : COST - (5*GEN[0,0] + 10*GEN[1,0] + 15*GEN[2,0] + 5*GEN[0,1] + 10*GEN[1,1] + 15*GEN[2,1] + 5*GEN[0,2] + 10*GEN[1,2] + 15*GEN[2,2] + 5*GEN[0,3] + 10*GEN[1,3] + 15*GEN[2,3] + 5*GEN[0,4] + 10*GEN[1,4] + 15*GEN[2,4] + 5*GEN[0,5] + 10*GEN[1,5] + 15*GEN[2,5]) : 1.0
demand_balance[0] : GEN[0,0] + GEN[1,0] + GEN[2,0] : 15.0
（略）

解へのアクセス２

添え字付き変数は、事後処理のために、解を配列やデータフレーム等に取り込むと便利です。変数$GEN$の解を配列sol[g,t]に入れるための書き方を比較します。

PuLP

value関数で変数を呼ぶと解を得ることができます。
変数との対応の辞書GEN(g,t)を作成していたので、それを利用して配列に値を代入します。

sol = np.array([[pulp.value(GEN[g, t]) for t in hours] for g in generators])

PYTHON-MIP

最もシンプルで、これだけで解の入った配列を得ることができます。

sol = GEN.astype(float) #ndarray

pyomo

記述方法はPuLPと似ています。なお、m1.GEN.getvalues()で、辞書形式の結果（キーが添え字の組み合わせ(g, t)、値が解）を得ることもできます。

sol = np.array([[pyo.value(m1.GEN[g, t]) for t in hours] for g in generators])

解のグラフ表示

必要な解を変数に代入すれば、後の処理はモデラーには関係なく行えます。

例えば配列sol[g,t]に入れた発電量の最適解をグラフ表示するには以下。

for g in generators:
  plt.bar(hours, sol[g], bottom=sol[:g].sum(axis=0), label='generator%d'%g)
plt.legend()
plt.xlabel('hour')
plt.ylabel('generation')
plt.show()

なお、上部で示した数値例に対する、発電量の最適解を図で示します。
最もコストが高いが発電量を柔軟に変化させられるgenerator2（緑）を活用するために、柔軟性の低いgenerator0（青）やgenerator1（橙）の発電量が抑えられていることがわかります。

おわりに

本記事では、"実用的"な変数・制約式から成る数理計画問題を対象として、Pythonで記述を行う方法を示しました。
dict型やリスト内包表記、演算子%など、Python特有の直観的にはわかりにくい記述方法もありますが、有料モデラー（GAMS等）を導入するまでもなく最適化をしてみたいという場合には、Pythonの各ライブラリを気軽に利用できます。