※記事内に商品プロモーションを含むことがあります。
はじめに
Pythonの最適化モデリングツールであるPyomoで、変数をベクトル化して扱う方法をまとめた。
変数をベクトル化することによって、変数の数が増えてもコード量の増加を抑えることができる。また、変数の数が不定な状況にも対応できる。
本記事では、0-1ナップサック問題を例としてコードを示す。
検証環境は以下の通り。
|
バージョン |
| Python |
3.7.4 |
| Pyomo |
5.6.9 |
| GLPK |
4.65 |
PyomoとGLPKのインストール方法は以下の記事を参照。
Pyomoで線形計画問題を解く
リンク
ナップサック問題
以下の0-1ナップサック問題を考える。
$$ \begin{array}{ll} \text{maximize} \ & \sum_{i=1}^{N} p_i x_i \\ \text{subject to} \ & \sum_{i=1}^{N} w_i x_i \le C \\ & x_i \in \{0, 1\}, (i=1, …, N) \end{array} $$
ここでは、$N=5$とし、
$\bm{p}=[2, 2, 5, 8, 12]$
$\bm{w}=[1, 2, 4, 6, 10]$
とする。
ナップサック問題の詳細は以下を参照
ナップサック問題と分枝限定法
Pyomoのソースコード
上記の線形計画問題をPyomoを使って記述し、最適化を実行したコードは以下のようになる。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
import pyomo.environ as pyo
prices = [2, 2, 5, 8, 12]
weights = [1, 2, 4, 6, 10]
capacity = 16
sense = pyo.maximize # pyo.minimize
N = len(prices)
N_index_list = list(range(N))
x_ini = [0]*N
xx = {i:x_ini[i] for i in range(N)}
def InitializeX(model, i):
"""変数の初期化"""
return xx[i]
def ObjRule(model):
"""目的関数"""
return sum(prices[i] * model.x[i] for i in N_index_list)
def ConstRule(model):
"""制約式"""
tmp = sum(weights[i] * model.x[i] for i in N_index_list)
if sense==pyo.maximize:
return tmp <= capacity
else:
return tmp >= capacity
model = pyo.ConcreteModel()
model.x = pyo.Var(N_index_list, domain=pyo.Binary,
initialize=InitializeX)
model.OBJ = pyo.Objective(rule=ObjRule, sense=sense)
model.Constraint = pyo.Constraint(rule=ConstRule)
opt = pyo.SolverFactory('glpk')
res = opt.solve(model, tee=False) # tee=Trueでソルバのメッセージを表示
print(f"評価関数：{model.OBJ()}")
for i in range(N):
print(f"x{i}: {model.x[i].value:.0f}")
print("選択された荷物")
print("price weight")
for i in range(N):
if model.x[i].value==1:
print(prices[i], weights[i])
|
実行結果は以下のようになる。得られた解は制約を満たしている。また、緩和問題は線形であるので、最適解（の1つ）が得られているはずである。
1
2
3
4
5
6
7
8
9
10
|
評価関数：20.0
x0: 0
x1: 0
x2: 0
x3: 1
x4: 1
選択された荷物
price weight
8 6
12 10
|
なお、コードをsense==pyo.minimizeと変更することで、最小化問題を解くことができる。このとき、制約式の符号は逆になる（≤から≥になる）。
ソースコードの解説
以下、ソースコードを簡単に解説する。
変数の定義
Varクラスを用いて変数を定義する。このとき、[0,1,...,N-1]というリストを引数に与えることで、ベクトル化された変数を定義できる。
また、関数InitializeXで変数の初期値を全て0で与えている（0-1ナップサック問題の場合、初期値を与えなくても収束への影響はほとんどない）。
1
2
3
4
5
6
7
8
9
10
11
12
|
N = len(prices)
N_index_list = list(range(N))
x_ini = [0]*N
xx = {i:x_ini[i] for i in range(N)}
def InitializeX(model, i):
"""変数の初期化"""
return xx[i]
...
model.x = pyo.Var(N_index_list, domain=pyo.Binary,
initialize=InitializeX)
|
目的関数の定義
次に、Objectiveクラスを用いて目的関数を定義する。ObjRule関数で目的関数を記述し、Objectiveクラスのruleに渡している。
なお、ObjRule関数の中身はジェネレータ式となっている（タプル内包表記とは呼ばない）。
1
2
3
4
5
|
def ObjRule(model):
"""目的関数"""
return sum(prices[i] * model.x[i] for i in N_index_list)
...
model.OBJ = pyo.Objective(rule=ObjRule, sense=sense)
|
制約条件の定義
Constraintクラスで制約条件を設定する。
1
2
3
4
5
6
7
8
9
|
def ConstRule(model):
"""制約式"""
tmp = sum(weights[i] * model.x[i] for i in N_index_list)
if sense==pyo.maximize:
return tmp <= capacity
else:
return tmp >= capacity
...
model.Constraint = pyo.Constraint(rule=ConstRule)
|
参考
Pyomoの公式リファレンス
Pyomo Documentation 5.7.2 — Pyomo 5.7.2 documentation
参考にさせていただいた日本語の記事
Pythonで線形計画問題を解く即戦力コード with Pyomo – Takala’s Memory
