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

はじめに

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)