はじめに #
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のコードは以下のようになります。
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}")
実行結果は以下のようになります。
評価関数：21.0
x0: 1
x1: 10
x2: 0
制約式の確認 #
制約式のdisplay()メソッドを使うと、最適化後の制約式の値や上限・下限を確認できます。
model.constraints.display()
実行結果
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通りの方法があります。
for i in range(N):
print(model.constraints[i].body())
# または
for c in model.constraints.values():
print(c.body())
制約式の上限値・下限値を直接取得したい場合、body()メソッドの代わりにub, lb属性にアクセスします。
実行結果（どちらも同じ）
12.0
10.0
2.0
ソースコードの解説 #
以下、制約式のベクトル化について、ソースコードを簡単に解説します。
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となります。
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)
