こんにちは！TECH Street編集部です。
TECH Streetでは、コミュニティメンバーが日々の学びや知見を「技術Tips」カテゴリで発信しています。
今回のテーマ
Pythonには無料の数理最適化ライブラリがあります。有名どころでは「OR-Tools」と「PuLP」が二大巨頭ではないでしょうか。
個人的にはOR-Tools推しですが、私の周囲ではPuLPを使用する人も多く、どちらが良いのか知りたくなりました。そこで、今回は「OR-Tools」と「PuLP」に数独を解かせて、その処理速度を比較して、どちらのライブラリがよいか検証してみます。
環境
問題設定
- 数独を用いる。
- 9x9の全マスに, 1～9の数字を必ず埋める。
- どの1行、どの1列、どの3x3においても同じ数字は1回だけしか使えない。
- 数字が埋まっている場所では、その数字を使う。
- 問題はこちらのページから借りてきました。
- 上級問題を10種類、各10回ずつ解かせ、時間を計測する。
定式化
- 定式化は整数計画法を用いた。
- 本来ならば目的関数を定義するが、今回の定式化では考慮しない。
- 今回の定式化はどちらかと言うとSATの考え方に近い。
- 9×9×9の箱を考える。
- 行方向×列方向×番号の3方向。
- 各セルに0か1の数値が入る。
- 変数は x(i, j, k)∈{0, 1}とする。
- 但し、i, j, k ∈ {0, 1, ..., 8}である。
- インデックス(i, j, k)はそれぞれ行方向、列方向、番号を表す。
- 例として x(1,1,7)=1なら、2行目、2列目のセルには8が入るという意味である。
- 既に数字が埋まっているインデックスと数値を以下で定義する。
- 例として、1行目、 1列目のセルに5が入っていたら、以下とする。
- 制約条件
- どの1行においても、同じ数字は1回だけしか使えない。
- どの1列においても、同じ数字は1回だけしか使えない。
- どの3x3においても、同じ数字は1回だけしか使えない。
- 数字が埋まっている場所では、その数字を使う。
実装
計算機のクラスを作成しています。重要なポイントは以下です。
- 計算機の抽象クラスには計算以外の機能を全て巻き取らせている。
- 時間計測
- 元になる配列のセッター、ゲッター
- 解のセッター、ゲッター
- PuLPの実装、 OR-Toolsの実装について
- 計算能力以外で差が出ないように、出来るだけ定式化のコードを同じにしている。
- それぞれの計算機で求解する場所の前後で時間計測を行っている。
結果
全体の結果は上図(箱ひげ図)の通りです。 平均で7倍前後程度の差が出ました。 また、計算一回当たりの処理時間の代表値は以下の表の通りです。
| 代表値 | OR-Tools | PuLP |
|---|---|---|
| 平均値[秒] | 0.013 | 0.099 |
| 中央値[秒] | 0.013 | 0.094 |
また、10回のトライアルの詳細な計算時間は上図の通りです。問題ごとに傾向はあまり見られない為、純粋なライブラリの実力が出ていると考えられます。
最後に
数独の求解に関して、7倍程度の差が出るとは予想外でした。
今回は変数が729個の整数計画問題として求解しましたので、 今後は線形計画法や混合整数計画法の求解でどの程度違うかを、 問題のスケールも含めて調べてみようかと思います。
Appendix1(定式化)
以下は定式化(数式)です。参考にしてください。
Appendix2(実装)
以下はクラスの全文です。 参考にしてください。
計算機の抽象クラスは以下です。
from abc import ABC, abstractmethod
from time import perf_counter_ns
from src.utils.enums.problem_types import ProblemTypes
class AbstractCalculator(ABC):
"""
計算機の抽象クラス
"""
@abstractmethod
def __init__(self) -> None:
"""
必要なら、イニシャライザは各サブクラスにて実装
"""
pass
@abstractmethod
def formulate(self) -> None:
"""
定式化処理は各サブクラスにて実装
"""
raise NotImplementedError("Implement each subclasses")
@abstractmethod
def calculate(self) -> None:
"""
求解処理は各サブクラスにて実装
"""
raise NotImplementedError("Implement each subclasses")
def set_array(self, array: list) -> None:
"""問題用の配列をセットする※初期化処理も兼ねる
Args:
array (list): 問題用の配列
"""
self.__array = array
self.__time_array = []
self.__calc_time = -1
self.set_answer_array([])
def set_answer_array(self, answer_array: list) -> None:
"""解答の配列をセットする
Args:
answer_array (list): 解答の配列
"""
self.__answer_array = answer_array
def set_problem_type(self, problem_type: ProblemTypes) -> None:
"""問題種別をセットする
Args:
problem_type (ProblemTypes): 問題種別
"""
self.__problem_type = problem_type
def append_time(self) -> None:
"""現在の時刻を配列に追加する
"""
self.__time_array.append(perf_counter_ns())
def set_calc_time(self) -> None:
"""時刻の配列から計算時間をセットする
Raises:
Exception: 時間計測していない
"""
time_array = self.__time_array
if len(time_array) < 2:
# 配列の長さが2より小さい->時間計測していない
raise Exception("No records")
else:
# 時刻の最初と最後を取得する
time0 = time_array[0]
time1 = time_array[1]
# ナノ秒で記録しているので, 秒単位に直す
calc_time = (time1-time0)*(10**-9)
# 0.1マイクロ秒より後の桁を切り捨てる
calc_time = round(calc_time, 7)
# 計算時間をセットする
self.__calc_time = calc_time
def get_array(self) -> list:
"""問題用の配列を出力する
Returns:
list: 問題用の配列
"""
return self.__array
def get_answer_array(self) -> list:
"""解答用の配列を出力する
Returns:
list: 解答用の配列
"""
return self.__answer_array
def get_problem_type(self) -> ProblemTypes:
"""問題種別を出力する
Returns:
ProblemTypes: 問題種別
"""
return self.__problem_type
def get_calc_time(self) -> float:
"""計算時間を出力する
Returns:
float: 計算時間
"""
return self.__calc_time
以下がPuLPの実装です。
import pulp
from src.calc.abstract_calculator import AbstractCalculator
from src.utils.enums.problem_types import ProblemTypes
class PulpCalculator(AbstractCalculator):
"""
PuLPを用いた計算機です。
"""
# region constructor
def __init__(self) -> None:
super().__init__()
# endregion constructor
# region formulate
def formulate(self) -> None:
"""
定式化処理。問題種別次第で処理を変更する。
"""
problem_type = self.get_problem_type()
if problem_type == ProblemTypes.IP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.LP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.SUDOKU:
self.__formulate_SUDOKU()
else:
raise ValueError(f"Invalid type: {problem_type}")
def __formulate_SUDOKU(self) -> None:
"""数独の定式化処理
"""
# モデルインスタンスを作成する
prob = pulp.LpProblem('Sudoku', pulp.LpMinimize)
# 目的関数はナシ
prob += 0
# 領域作成
self.__numbers = list(range(1, 10))
self.__xs = list(range(1, 10))
self.__ys = list(range(1, 10))
numbers = self.__numbers
xs = self.__xs
ys = self.__ys
choices = pulp.LpVariable.dicts(
"Choice",
(xs, ys, numbers),
0,
1,
pulp.LpInteger
)
# 制約(k方向にどれか一個1)(答えの1マスには一つの数値しか存在できない）
for i in xs:
for j in ys:
prob += (sum(choices[i][j][k] for k in numbers) == 1)
# 制約(i方向にどれか一個1)(答えの1行には各数値は一つずつしか存在できない）
for j in ys:
for k in numbers:
prob += (sum(choices[i][j][k] for i in xs) == 1)
# 制約(j方向にどれか一個1)(答えの1列には各数値は一つずつしか存在できない）
for i in xs:
for k in numbers:
prob += (sum(choices[i][j][k] for j in ys) == 1)
# 制約(3x3に区切られたマス内でも各数字は一つしか存在できない)
for k in numbers:
for i3 in range(3):
prob += (
sum(
choices[i3*3+i][j][k]
for i in range(1, 4) for j in range(1, 4)
) == 1
)
prob += (
sum(
choices[i3*3+i][j][k]
for i in range(1, 4) for j in range(4, 7)
) == 1
)
prob += (
sum(
choices[i3*3+i][j][k]
for i in range(1, 4) for j in range(7, 10)
) == 1
)
# 制約(既に数値が置いてある場所に変数をセットしていく)
array = self.get_array()
for i, row in enumerate(array):
for j, cell in enumerate(row):
if cell == 0:
# 数値が確定していないので次へ
pass
else:
# 数値が1~9のどれかが入っている
prob += (choices[i+1][j+1][cell] == 1)
self.__choices = choices
self.__prob = prob
# endregion formulate
# region calculate
def calculate(self) -> None:
"""計算機処理. 問題種別次第で処理を変更する。
"""
problem_type = self.get_problem_type()
if problem_type == ProblemTypes.IP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.LP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.SUDOKU:
self.__calculate_SUDOKU()
else:
raise ValueError(f"Invalid type: {problem_type}")
def __calculate_SUDOKU(self):
"""数独の求解処理
"""
# モデルインスタンスを取得する
prob = self.__prob
self.append_time()
# 求解処理を行うが, メッセージが鬱陶しいので見せないことにする
state = prob.solve(pulp.PULP_CBC_CMD(msg=False))
self.append_time()
self.set_calc_time()
# 結果を解答用配列にセットする
if pulp.LpStatus[state] == "Optimal":
numbers = self.__numbers
xs = self.__xs
ys = self.__ys
choices = self.__choices
answer_array = []
for x in xs:
row = []
for y in ys:
for v in numbers:
if choices[x][y][v].value() == 1:
row.append(v)
answer_array.append(row)
self.set_answer_array(answer_array)
else:
print("Couldn't optimize")
# endregion calculate
以下がOR-Toolsの実装
from ortools.linear_solver import pywraplp
from src.calc.abstract_calculator import AbstractCalculator
from src.utils.enums.problem_types import ProblemTypes
class OrToolsCalculator(AbstractCalculator):
"""
OR-Toolsを用いた計算機
"""
# region constructor
def __init__(self) -> None:
super().__init__()
# endregion constructor
# region formulate
def formulate(self) -> None:
"""
定式化処理。問題種別次第で処理を変更する。
"""
problem_type = self.get_problem_type()
if problem_type == ProblemTypes.IP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.LP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.SUDOKU:
self.__formulate_SUDOKU()
else:
raise ValueError(f"Invalid type: {problem_type}")
def __formulate_SUDOKU(self) -> None:
"""数独の定式化処理
"""
# モデルインスタンスを作成する
solver: pywraplp.Solver = pywraplp.Solver.CreateSolver("SAT")
# 領域作成
var_dict = {}
for i in range(9): # 行方向
for j in range(9): # 列方向
for k in range(9): # 奥行方向
var_name = f"ans_{i}_{j}_{k}"
var_dict[(i, j, k)] = solver.IntVar(0, 1, var_name)
# 制約(k方向にどれか一個1)(答えの1マスには一つの数値しか存在できない）
for i in range(9):
for j in range(9):
solver.Add(sum(var_dict[(i, j, k)] for k in range(9)) == 1)
# 制約(i方向にどれか一個1)(答えの1行には各数値は一つずつしか存在できない）
for j in range(9):
for k in range(9):
solver.Add(sum(var_dict[(i, j, k)] for i in range(9)) == 1)
# 制約(j方向にどれか一個1)(答えの1列には各数値は一つずつしか存在できない）
for i in range(9):
for k in range(9):
solver.Add(sum(var_dict[(i, j, k)] for j in range(9)) == 1)
# 制約(3x3に区切られたマス内でも各数字は一つしか存在できない)
for k in range(9):
for i3 in range(3):
solver.Add(
sum(
var_dict[(i3*3+i, j, k)]
for i in range(3) for j in range(0, 3)
) == 1
)
solver.Add(
sum(
var_dict[(i3*3+i, j, k)]
for i in range(3) for j in range(3, 6)
) == 1
)
solver.Add(
sum(
var_dict[(i3*3+i, j, k)]
for i in range(3) for j in range(6, 9)
) == 1
)
# 制約(既に数値が置いてある場所に変数をセットしていく)
array = self.get_array()
for i, row in enumerate(array):
for j, cell in enumerate(row):
if cell == 0:
# 数値が確定していないので次へ
pass
else:
# 数値が1~9のどれかが入っている
# 1を引いてインデックスを取る
k = cell - 1
solver.Add(var_dict[(i, j, k)] == 1)
# 変数のディクショナリとモデルインスタンスをセットする
self.__var_dict = var_dict
self.__solver = solver
# endregion formulate
# region calculate
def calculate(self) -> None:
"""計算機処理。問題種別次第で処理を変更する。
"""
problem_type = self.get_problem_type()
if problem_type == ProblemTypes.IP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.LP:
raise NotImplementedError("Not implemented")
elif problem_type == ProblemTypes.SUDOKU:
self.__calculate_SUDOKU()
else:
raise ValueError(f"Invalid type: {problem_type}")
def __calculate_SUDOKU(self):
"""数独の求解処理
"""
# ソルバを取得する
solver: pywraplp.Solver = self.__solver
# 求解する
self.append_time()
status = solver.Solve()
self.append_time()
self.set_calc_time()
# 結果をセットする
if ((status == pywraplp.Solver.OPTIMAL) or (status == pywraplp.Solver.FEASIBLE)):
# 解答用の配列のひな形を作成する
answer_array = [
[0 for _ in range(9)]
for _ in range(9)
]
var_dict = self.__var_dict
for i in range(9):
for j in range(9):
for k in range(9):
# 変数インスタンスを取り出し、インデックス+1をして解答用の配列にセットする
variable: pywraplp.Variable = var_dict[(i, j, k)]
val = variable.solution_value()
if val != 0:
answer_array[i][j] = k+1
break
# 解答用の配列をセットする
self.set_answer_array(answer_array)
else:
print("Couldn't optimize")
# endregion calculate
記事執筆者
小澤 陽介さん
パーソルキャリア株式会社
テクノロジー本部 デジタルテクノロジー統括部 デジタルビジネス部 アナリティクスグループ
前々職は半導体製造会社にて、リソグラフィ工程のプロセスエンジニアリングに従事。前職では機械系メーカーR&Dにて数理最適化のアプリケーション開発及び工程可視化に従事。パーソルキャリアに入ってからはアナリストとして、主に生成系AIに関連するアプリケーション開発や検証に従事。意思決定の為の計算機開発にも携わっている。
