数理最適化モデルの要素をクラスとして実装するプラクティス【Pyomo編】

※記事内に商品プロモーションを含むことがあります。

はじめに

複数の要素から構成される最適化モデルを最適化ライブラリPyomoで実装するとき、コードを再利用しやすくするため、各要素をクラスとして実装する方法を検討しました。

複数の要素から構成されるモデルとは以下のようなものです。

ビルの空調システム（個々の部屋や空調装置が要素）
発送電システム（個々の発電装置や受電装置、蓄電装置が要素）

要素をクラスとして実装することにより、異なる構成のシステムにコードを再利用しやすくなります。

検討に使用したPythonとPyomoのバージョンは以下の通りです。

ソフトウェア	バージョン
Python	3.9.7
Pyomo	6.4.1

CasADi+Python版の記事は以下です。
数理最適化モデルの要素をクラスとして実装するプラクティス【CasADi+Python編】 – Helve Tech Blog

例題とするシステム

例題とするシステムは、画像のように4つのプロセス（要素）から構成されています。各プロセスを通過する量$x_1, …, x_4$が大きいほど利益が多くなるとします。また、$x_1, …, x_4$にはそれぞれ上限があるとします。

process-network

最適化問題を次式のように定義します。システムを通過する量の和は6とします（最後の制約式）。

$$ \max \ x_1 + 2x_2 + 3x_3 + x_4 \\ \mathrm{s.t.} \ x_1 \le 5 \\ x_2 \le 3 \\ x_3 \le 3 \\ x_4 \le 5 \\ x_1 + x_2 = x_3 + x_4 \\ x_1 + x_2 = 6 $$

要素をクラス化した実装例

プロセスをProcessクラスとして実装した最適化のコードを以下に示します。__init__メソッドで変数や制約式を定義します。

proc1.xやproc2.profitのように、最適化の変数をインスタンス変数として扱うことができます。

 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


import pyomo.environ as pe

class Process:
    count = 0
    def __init__(self, mdl, x_ub, unit_profit):
        # オブジェクト名、変数、パラメータ、制約式を設定する
        name = "proc"+str(Process.count)
        self.name = name
        Process.count += 1

        exec(f"mdl.{name}_x = pe.Var(domain=pyo.Reals, bounds=(0, x_ub))")
        exec(f"self.x = mdl.{name}_x")

        exec(f"mdl.{name}_unit_profit = pe.Param(default=unit_profit)")
        exec(f"self.unit_profit = mdl.{name}_unit_profit")

        exec(f"mdl.{name}_profit = pe.Expression(expr = mdl.{name}_unit_profit*mdl.{name}_x)")
        exec(f"self.profit = mdl.{name}_profit")

mdl = pe.ConcreteModel()
proc1 = Process(mdl, 5, 1)
proc2 = Process(mdl, 3, 2)
proc3 = Process(mdl, 3, 3)
proc4 = Process(mdl, 5, 1)

mdl.const1 = pe.Constraint(expr = proc1.x + proc2.x == proc3.x + proc4.x)
mdl.const2 = pe.Constraint(expr = proc1.x + proc2.x == 6)

mdl.obj = pe.Objective(expr=proc1.profit + proc2.profit + proc3.profit + proc4.profit,
                       sense = pe.maximize)

opt = pe.SolverFactory('glpk')
res = opt.solve(mdl, tee=False)

print(f"Objective: {mdl.obj()}")

for v in mdl.component_data_objects(ctype=pe.Var):
    print(f"{v.name=}, {v.lb=}, {v.value=}, {v.ub=}")

実行結果

実行すると、以下のように最適化結果が表示されます。

1
2
3
4
5


Objective: 21.0
v.name='proc1_x', v.lb=0, v.value=3.0, v.ub=5
v.name='proc2_x', v.lb=0, v.value=3.0, v.ub=3
v.name='proc3_x', v.lb=0, v.value=3.0, v.ub=3
v.name='proc4_x', v.lb=0, v.value=3.0, v.ub=5

例題では4つのプロセスの構造（制約式やパラメータ）が同じであるため、クラスを1つしか作りませんが、異なる構造のプロセスがある場合、その分だけ対応したクラスを作成します。

要素クラスの解説

要素のProcessクラスを解説します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


class Process:
    count = 0
    def __init__(self, mdl, x_ub, unit_profit):
        # 変数、パラメータ、制約式を設定する
        name = "proc"+str(Process.count)
        self.name = name
        Process.count += 1

        exec(f"mdl.{name}_x = pe.Var(domain=pe.Reals, bounds=(0, x_ub))")
        exec(f"self.x = mdl.{name}_x")

        exec(f"mdl.{name}_unit_profit = pe.Param(default=unit_profit)")
        exec(f"self.unit_profit = mdl.{name}_unit_profit")

        exec(f"mdl.{name}_profit = pe.Expression(expr = mdl.{name}_unit_profit*mdl.{name}_x)")
        exec(f"self.profit = mdl.{name}_profit")

変数self.nameはオブジェクトの名前（プロセスの固有の名前）です。self.nameはオブジェクト固有の変数名を付けるために使用します。

Pyomoでは、最適化の変数やパラメータ、制約式をModelオブジェクト（コードでは変数mdl）の変数として定義します。すなわち、同じクラスで複数の要素（プロセス）を定義する場合、変数名が重複しないようにします。そのため、exec関数を使っています。

Varクラスで最適化の変数を定義しています。Realsは実数であることを設定し、boundsは変数の上限・下限で指定します。

Paramでパラメータ（最適化計算中も一定値の変数）を定義します。default引数に値を渡します。

Expressionクラスで式を定義します。ここでは、{name}_profit = mdl.{name}_unit_profit*mdl.{name}_xという関係式を定義しています。

また、mdlの変数をオブジェクト変数にも格納しています。これは、評価関数などを定義するときに、オブジェクトの外側から変数にアクセスしやすくするためです。

クラス化しない場合のコード

クラス化しない場合のコード例を示します。変数をベクトル化した場合、もう少し簡潔に記述できます。

 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


mdl = pe.ConcreteModel()

mdl.proc1_x = pe.Var(domain=pe.Reals, bounds=(0, 5))
mdl.proc1_unit_profit = pe.Param(default=1)
mdl.proc1_profit = pe.Expression(expr = mdl.proc1_unit_profit*mdl.proc1_x)

mdl.proc2_x = pe.Var(domain=pe.Reals, bounds=(0, 3))
mdl.proc2_unit_profit = pe.Param(default=2)
mdl.proc2_profit = pe.Expression(expr = mdl.proc2_unit_profit*mdl.proc2_x)

mdl.proc3_x = pe.Var(domain=pe.Reals, bounds=(0, 3))
mdl.proc3_unit_profit = pe.Param(default=3)
mdl.proc3_profit = pe.Expression(expr = mdl.proc3_unit_profit*mdl.proc3_x)

mdl.proc4_x = pe.Var(domain=pe.Reals, bounds=(0, 5))
mdl.proc4_unit_profit = pe.Param(default=1)
mdl.proc4_profit = pe.Expression(expr = mdl.proc4_unit_profit*mdl.proc4_x)

mdl.const1 = pe.Constraint(expr = mdl.proc1_x + mdl.proc2_x\
                           == mdl.proc3_x + mdl.proc4_x)
mdl.const2 = pe.Constraint(expr = mdl.proc1_x + mdl.proc2_x == 6)

mdl.obj = pe.Objective(expr=mdl.proc1_profit+mdl.proc2_profit+mdl.proc3_profit+mdl.proc4_profit,
                       sense = pe.maximize)

opt = pe.SolverFactory('glpk')
res = opt.solve(mdl, tee=False)

print(f"Objective: {mdl.obj()}")

for v in mdl.component_data_objects(ctype=pe.Var):
    print(f"{v.name=}, {v.lb=}, {v.value=}, {v.ub=}")