Pyomoの制約をベクトル化する

$$ \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