この記事は数理最適化Advent Calender 2022の10日目の記事です。

1.はじめに

書籍「Pythonではじめる数理最適化」¹の著者の岩永です。出版して1年が過ぎましたが、多くの読者と支援者に恵まれ、先日第5刷を発行できました。ご支援くださったみなさまに感謝申し上げます。

本記事では、書籍「Pythonではじめる数理最適化」をテーマに開催したORセミナー²の参加者アンケートで「勉強になった！」と評判が良かったクラス設計について取り上げます。本記事を読むことで数理最適化プログラムのクラス設計の型とポイントについて理解することができます。また、モデリング言語はPythonライブラリのPuLPを利用します。予めパッケージ管理ツールpipやcondaなどを用いてインストールしておいてください。

2.クラス化は処理のメソッド化から

はじめにお断りしておくと、クラスは使いますがオブジェクト指向は採用しません。数理最適化プログラムは、本質的に関数を直列に繋ぐタイプの処理で、順序依存が強いためオブジェクト指向とは相性が悪いためです。しかしながら、処理の過程で常に状態を把握したいといった要望や、保守運用性を高めたいといった要望があるためクラスの機能を使います³。
数理最適化プログラムをクラス化する１番のメリットは、各処理をメソッドに切り分けることで、適度な分量のコードに分割でき、可読性が高まることです。例えば、数理最適化プログラムの決まった処理に「・・・→集合・定数の作成→モデリング→最適化計算→・・・」といった一連の流れがあります。クラス化のはじめの一歩はこの処理をメソッドに分割することです。

3.数理最適化モデル

本記事では、生産計画問題の数理最適化モデルを扱いますが、クラス設計が主題なのでコードを読みながら必要に応じて参照してかまいません。書籍では第2章の内容になります。

集合
- 製品の集合：$P$
- 原料の集合：$M$
定数
- 在庫量：$stock_m\ \ (m\in M)$
- 必要量：$require_{p,m}\ \ (p\in P, m\in M)$
- 利得：$gain_p\ \ (p\in P)$
変数
- 生産量：$x_p \in \mathbb{R}\ \ (p\in P)$
制約式
- 製品の生産量は$0$以上：
  $x_p \geq 0\ \ (\forall p\in P)$
- 生産は在庫の範囲で行う：
  $\sum_{p\in P}require_{p,m}\cdot x_{p,m} \leq stock_m\ \ (\forall m\in M)$
目的関数(最大化)
利得の合計：$\sum_{p \in P}gain_p \cdot x_p$

4.コード

さっそくですがコードの全容を掲載します。50行程度の小さなプログラムなので上から読んでいっても読み切れますが、通常はif __name__ == '__main__':のブロックから読み始めます。コードの後で解説を行います。

import pulp

class ProdPlan:
	def __init__(self, P, M, m2s, p2g, pm2r):
        # 集合と定数
		self.P = P
		self.M = M
		self.m2s = m2s
		self.p2g = p2g
		self.pm2r = pm2r
        
        # 数理最適化モデル
		self.x = None
		self.model = None

        # 最適化計算ステータス
		self.status = None
		return

	def modeling(self):
		# 数理モデル作成
		self.model = pulp.LpProblem('ProdPlan', pulp.LpMaximize)

		# 変数
		self.x = pulp.LpVariable.dicts('x', self.P, cat='Continuous', lowBound=0)

		# 制約式
		for m in self.M:
			self.model += pulp.lpSum(self.pm2r[p,m] * self.x[p] for p in self.P) <= self.m2s[m]

		# 目的関数
		self.model += pulp.lpSum(self.p2g[p] * self.x[p] for p in self.P)
		return

	def solve(self):
		self.status = self.model.solve()
		return

if __name__ == '__main__':
	# (1)集合・定数の作成
	P = ['p1', 'p2', 'p3', 'p4']
	M = ['m1', 'm2', 'm3']
	m2s = {'m1': 35, 'm2': 22, 'm3': 27}
	p2g = {'p1': 3, 'p2': 4, 'p3': 4, 'p4': 5}
	pm2r = {('p1', 'm1'): 2, ('p1', 'm2'): 0, ('p1', 'm3'): 1
			, ('p2', 'm1'): 3, ('p2', 'm2'): 2, ('p2', 'm3'): 0
			, ('p3', 'm1'): 0, ('p3', 'm2'): 2, ('p3', 'm3'): 2
			, ('p4', 'm1'): 2, ('p4', 'm2'): 2, ('p4', 'm3'): 2}

	# (2-1)モデリング（集合・定数の設定）
	prod_plan = ProdPlan(P, M, m2s, p2g, pm2r)

	# (2-2)モデリング（変数・制約・目的関数の設定）
	prod_plan.modeling()

	# (3)最適化計算
	prod_plan.solve()

5.コードの解説

メインルーチン

まず、if __name__ == '__main__':のブロックから確認します。

if __name__ == '__main__':
	# (1)集合・定数の作成
	P = ['p1', 'p2', 'p3', 'p4']
	M = ['m1', 'm2', 'm3']
	m2s = {'m1': 35, 'm2': 22, 'm3': 27}
	p2g = {'p1': 3, 'p2': 4, 'p3': 4, 'p4': 5}
	pm2r = {('p1', 'm1'): 2, ('p1', 'm2'): 0, ('p1', 'm3'): 1
			, ('p2', 'm1'): 3, ('p2', 'm2'): 2, ('p2', 'm3'): 0
			, ('p3', 'm1'): 0, ('p3', 'm2'): 2, ('p3', 'm3'): 2
			, ('p4', 'm1'): 2, ('p4', 'm2'): 2, ('p4', 'm3'): 2}

	# (2-1)モデリング（集合・定数の設定）
	prod_plan = ProdPlan(P, M, m2s, p2g, pm2r)

	# (2-2)モデリング（変数・制約・目的関数の設定）
	prod_plan.modeling()

	# (3)最適化計算
	prod_plan.solve()

(1)(2-1)(2-2)(3)のコメントから、それぞれの処理が数理最適化問題を解く一般的な流れになっていることを確認できます。

(1)集合・定数の作成
(2-1)モデリング（集合・定数の設定）
(2-2)モデリング（変数・制約・目的関数の設定）
(3)最適化計算

このように処理をメソッドに分けるだけで全体の見通しがよくなります。さらにこのようにメソッドに分けることでデバッグしやすくなることが重要です。数理最適化プログラムの実装では、通常のプログラム開発と比べてバグが頻発します。その原因は次の3点に代表されます。

データに誤りがある
数理モデルに誤りがある
最適化計算に誤りがある

データ作成、モデリング、最適化計算の処理が分離されていないとデバッグが難しく、調査に何倍も時間がかかります。上記のように(1)(2-1)(2-2)(3)のように分離することでデバッグプリントを挟めるのも利点です。例えば、処理の前後に数理最適化モデルの状態を出力する次のメソッドを実装しておくことをおすすめします。

集合・定数を設定した後に、データを表示するshow_dataメソッドを入れる
モデリングした後に、数理最適化モデルを表示するshow_modelメソッドを入れる
最適化計算した後に、最適化結果を表示するshow_resultメソッドを入れる

(1)集合・定数の作成

メインルーチンの次のコードで集合・定数の作成を行っています。

	P = ['p1', 'p2', 'p3', 'p4']
	M = ['m1', 'm2', 'm3']
	m2s = {'m1': 35, 'm2': 22, 'm3': 27}
	p2g = {'p1': 3, 'p2': 4, 'p3': 4, 'p4': 5}
	pm2r = {('p1', 'm1'): 2, ('p1', 'm2'): 0, ('p1', 'm3'): 1
			, ('p2', 'm1'): 3, ('p2', 'm2'): 2, ('p2', 'm3'): 0
			, ('p3', 'm1'): 0, ('p3', 'm2'): 2, ('p3', 'm3'): 2
			, ('p4', 'm1'): 2, ('p4', 'm2'): 2, ('p4', 'm3'): 2}

このコードのポイントは数理モデルの集合と定数はPythonではリストと辞書を使って定義する点です。数理モデルの集合をPythonでも集合（set）を使って実装したくなりますが、一般的に集合は順序が不定であるため弊害があります。例えば、次の記事「整数計画ソルバーが実行毎に異なる最適解を出力する謎を解け！」で問題点が指摘されています。数学上で同値だとしても、ソルバーの実装上では同値にならないことがあります。

また、上記のプログラムではリストと辞書の内容をベタ書きしましたが、データを取得する際にデータ解析ライブラリpandasのread_csv関数を利用するケースが多くあります。pandasを用いてデータ取得した場合もデータをリストと辞書に変換しておくと便利です。pandasのデータフレームを利用して、(2-2)のモデリングをすることもできますが、次の2点の理由でリストと辞書に変換しておくことをおすすめします。
・PuLPによる数理モデリングはPythonのリストと辞書を利用すると自然な記述ができ、可読性が高い
・pandasなどのライブラリに依存せず、組み込まれている機能のみで記述できる

(2-1)モデリング（集合・定数の設定）

メインルーチンのprod_plan = ProdPlan(P, M, m2s, p2g, pm2r)で集合と定数の設定をしています。def __init__(self, P, M, m2s, p2g, pm2r):のブロックを確認します。

	def __init__(self, P, M, m2s, p2g, pm2r):
        # 集合と定数
		self.P = P
		self.M = M
		self.m2s = m2s
		self.p2g = p2g
		self.pm2r = pm2r
        
        # 数理最適化モデル
		self.x = None
		self.model = None

        # 最適化計算ステータス
		self.status = None
		return

ここでは、数理最適化モデルの集合と辞書を設定しています。また、数理最適化モデルや最適化計算後の情報も定義されていることに注意してください。個人的には次の2つの情報を含める場合があります。

最適化計算のオプション(PuLPのオプションであるtimeLimit、gapRelなど。)
最適化計算の結果(目的関数値、変数値から計算される各種値)

(2-2)モデリング（変数・制約・目的関数の設定）

メインルーチンのprod_plan.modeling()で数理最適化モデルを定義します。クラスのメソッドであるdef modeling(self):のブロックを確認します。


	def modeling(self):
		# 数理モデル作成
		self.model = pulp.LpProblem('ProdPlan', pulp.LpMaximize)

		# 変数
		self.x = pulp.LpVariable.dicts('x', self.P, cat='Continuous', lowBound=0)

		# 制約式
		for m in self.M:
			self.model += pulp.lpSum(self.pm2r[p,m] * self.x[p] for p in self.P) <= self.m2s[m]

		# 目的関数
		self.model += pulp.lpSum(self.p2g[p] * self.x[p] for p in self.P)
		return

数理モデルと変数を作成する箇所self.model=...、self.x=...で数理最適化モデルの情報を保存していることに注意してください。また、変数は、定数と同様にPythonの辞書で定義しています。そうすることで制約式と目的関数が自然に定義でき、可読性が高いことを確認できます。例えば、

for m in self.M:
	self.model += pulp.lpSum(self.pm2r[p,m] * self.x[p] for p in self.P) <= self.m2s[m]

というコードは、

各原料について、原料の利用量の総和は在庫量以下
ただし、原料の利用量は、製品の生産量と原料の必要量を掛け合わせたもの

と自然に解釈することができます⁴。モデリング言語を用いる場合、できるだけ行列やベクトル表現を避けることをおすすめします。行列やベクトルは便利で簡潔な記述を提供しますが、式の意味が隠蔽されるため可読性が低かったり、拡張性や柔軟性も損いやすいためです⁵。行列やベクトルを使った数理最適化モデルは、自分で実装したものはわかりやすいのですが、他人の実装を読むと「ちょっとわかりにくい」と思うことがちらほらあるのではないでしょうか。また、仕様変更により制約の追加や変更をすると「あ、面倒くさい」と気づくはずです。モデリング言語は、より多くの人が効率よく数理最適化の技術を利用できることを目的に開発されたものなので、可読性の高さと拡張性・柔軟性を担保することは理にかなっています。

(3)最適化計算

メインルーチンのprod_plan.solve()で最適化計算を実行します。クラスのメソッドであるdef solve(self):のブロックを確認します。

	def solve(self):
		self.status = self.model.solve()

ここでは、最適化計算を実行するself.status = self.model.solve()の1行しかありませんが、数理最適化プログラムをモジュール化していくと最適化実行部分を分離しておくと便利な場合があります。個人的には次の3つの情報を含める場合があります。

ソルバーの選択
最適化オプションの指定
最適化計算結果から計算される各種値を数理モデルに保存

また、上記のコードでは実行結果の確認をしていませんが、もし最適化結果を表示させたい場合、次のshow_resultメソッドなどを実装してメインルーチンで確認するとよいでしょう。

    def show_result(self):
		print('Status:', pulp.LpStatus[self.status])
		for p in self.P:
			print(p, self.x[p].value())

6.おわりに

本記事では、数理最適化プログラムのクラス設計の一例を紹介しました。本来、クラス設計は自由であるものですが、数理最適化プログラムの特徴を考慮するとテンプレートがあってもよいと思います。
数理最適化プログラムはバグが頻発し、その原因の調査が大変であることを記事で述べました。データの誤り、数理最適化モデルの誤り、最適化計算の誤りなどがありますが、プログラム的な難しさだけでなく、数学的な難しさにも向き合わないと原因究明ができません。そのため、クラス設計はデバッグのしやすさを重要視すべきです。
また、数理最適化プログラムを社会実装するシーンでは仕様変更が頻繁に起こるので、プログラムの拡張性や柔軟性を担保することも重要です。クラスの機能をフルに使いこなせるようになれば、可読性が向上するだけでなく、プログラムの汎用性も高まります。すると、エラーハンドリングやテストの実装、実験の自動化など様々な利点が生まれます。

本記事が数理最適化の初学者にとって有益な情報になれば幸いです。もし、「Pythonではじめる数理最適化」の続本が出せるなら本記事の内容も書きますね。いつになることやら、、、

リンク集

「Pythonではじめる数理最適化」のオーム社のページ
「Pythonではじめる数理最適化」のAmazonのページ
「Pythonではじめる数理最適化の全コードが公開されているGitHub
ORセミナー
オペレーションズ・リサーチ学会
モデリング言語：PuLP

「Pythonではじめる数理最適化」は、岩永二郎、石原響太、西村直樹、田中一樹の4人の著者により執筆し、2021年9月21日にオーム社から出版されています。 ↩
ORセミナーは、オペレーションズ・リサーチ学会主催のセミナーです。 ↩
2022年12月13日追記部分。数理最適化プログラムをオブジェクト指向で実装した場合、どこかで破綻するか、保守運用したくなくなるくらい複雑なプログラムになります。また、本記事で紹介する方法も完全ではありません。大規模数理最適化問題では、計算やメモリの効率化の問題でデータ作成部分とモデリング部分を切り離せなくなる場合もあり、結局クラスの設計が破綻することがあります。 ↩
わかりやすさを重視して省略した解釈をしています。 ↩
数理最適化の理論やアルゴリズムを学ぶためには行列やベクトルを用いた表現に慣れておく必要があります。数理最適化問題を表現する立場と数理最適化問題を解く立場では何が自然で何が便利なのかが異なります。 ↩