※記事内に商品プロモーションを含むことがあります。
はじめに
Pythonの最適化モデリングツールであるPyomoで、制約をベクトル化して扱う方法をまとめました。変数をベクトル化する記事は以下です。
Pyomoの変数をベクトル化する – Helve Tech Blog
制約をベクトル化することによって、制約の数が増えてもコード量の増加を抑えることができます。また、制約の数が不定な状況にも対応できます。
検証環境は以下の通りです。GLPKという線形ソルバを使用します。
- Python 3.9.7
- Pyomo 6.4.1
- GLPK 4.65
PyomoとGLPKのインストール方法は以下の記事を参照ください。
Pyomoで線形計画問題を解く – Helve Tech Blog
リンク
線形計画問題
以下の線形計画問題を考えます。
$$ \begin{array}{ll} \text{maximize} \ & x_0 + 2x_1 + 3x_2 \\ \text{subject to} \ & 2x_0 + x_1 \le 12 \\ & x_1 + 3x_2 \le 10 \\ & x_2 + 2x_0 \le 8 \\ & x_0, x_1, x_2 > 0 \end{array} $$
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
|
import pyomo.environ as pyo
N = 3 # 変数の数
obj_coef = [1, 2, 3]
coef1 = [2, 1, 1]
coef2 = [1, 3, 2]
intercept = [12, 10, 8]
def obj_rule(model):
"""目的関数"""
return sum(obj_coef[i]*model.x[i] for i in range(N))
def const_rule(model, i):
"""制約式"""
return coef1[i]*model.x[i]+coef2[i]*model.x[(i+1)%N] <= intercept[i]
model = pyo.ConcreteModel()
model.x = pyo.Var(range(N), domain=pyo.NonNegativeReals)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)
model.constraints = pyo.Constraint(range(N), rule=const_rule)
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}")
|
実行結果は以下のようになります。
1
2
3
4
|
評価関数：21.0
x0: 1
x1: 10
x2: 0
|
制約式の確認
制約式のdisplay()メソッドを使うと、最適化後の制約式の値や上限・下限を確認できます。
1
|
model.constraints.display()
|
実行結果
1
2
3
4
5
|
constraints : Size=3
Key : Lower : Body : Upper
0 : None : 12.0 : 12.0
1 : None : 10.0 : 10.0
2 : None : 2.0 : 8.0
|
また、最適化後の制約式の値を直接取得するには、2通りの方法があります。
1
2
3
4
5
|
for i in range(N):
print(model.constraints[i].body())
# または
for c in model.constraints.values():
print(c.body())
|
制約式の上限値・下限値を直接取得したい場合、body()メソッドの代わりにub, lb属性にアクセスします。
実行結果（どちらも同じ）
ソースコードの解説
以下、制約式のベクトル化について、ソースコードを簡単に解説します。
1
2
3
4
5
|
def const_rule(model, i):
"""制約式"""
return coef1[i]*model.x[i] + coef2[i]*model.x[(i+1)%N] <= intercept[i]
model.constraints = pyo.Constraint(range(N), rule=const_rule)
|
制約式をconst_ruleという関数で定義しています。最初の引数はPyomoの最適化モデルのオブジェクト、2番目の引数は制約式のインデックスです。(i+1)%Nは少しトリッキーですが、i=2のときインデックスを0に戻すための工夫です。
次に、Constraintクラスで複数の制約式を同時に定義します。このクラスの最初の引数は制約式のインデックスです。制約式が3つあるのでrange(N)としています。また、rule引数に制約式を定義する関数を与えます。
より複雑な制約の与え方
詳細は触れませんが、以下のように制約式のインデックスを増やすことで、より複雑な制約式を定義できます。このとき、制約式の数はN*Mとなります。
1
2
3
4
5
|
def const_rule2(model, i, j):
"""制約式"""
return coef1[i]*model.x[j] + coef2[i]*model.x[(j+1)%N] <= intercept[i]
model.constraints2 = pyo.Constraint(range(N), range(M), rule=const_rule2)
|
参考
Pyomoの公式リファレンス
Constraints — Pyomo 6.4.1 documentation
