JavaScriptを有効にしてください

Pyomoの変数をベクトル化する

 ·  ☕ 4 min read  ·  ✍️ Helve

はじめに

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

シェアする
ブログランキングに参加中です。クリックして頂けると励みになります。

Helve
WRITTEN BY
Helve
関西在住、電機メーカ勤務のエンジニア。Twitterで新着記事を配信中です