遺伝的アルゴリズムで難しいインスタンスを生成してみた
この記事は数理最適化 Advent Calendar 2024の21日目の記事です。
導入
数理最適化の研究や実装において、作成したアルゴリズムの性能検証や異なる手法間のベンチマークを取得するのに使うためのテスト用のインスタンスが欲しい状況が往々にして生じる。
しかしながら、乱数を用いてインスタンスを生成するとしばしば簡単すぎる問題が生成されがちであるため、「そこそこ難しい」 インスタンスを用意するのは意外と厄介なタスクである。
ベンチマークを取得する対象に任意性があるのであればMIPLIB等のインスタンスを用いてベンチマークを取れば良いものの、自作した数理モデルでベンチマークを取る場合には上記の困難は避けられない。
本記事では、遺伝的アルゴリズムを用いて恣意的に難しいインスタンスを生成してみようと思う。
また、今回のコード全てはNotebookとしてGitHub上に公開しており、BinderやColab環境で試すことが可能である。
インスタンス生成戦略
インスタンスの「難しさ」の定義にはいろいろあるが、一番わかりやすい定義の一つとして「ソルバーの求解時間」が挙げられる。一般に計算時間が長ければ長いほどソルバーが上下界のギャップを埋めるのに苦労しているといえるので、最適解を求めるのが「難しい」インスタンスと言えるだろう。[1]
そのため、
- まず乱数を用いてランダムにインスタンスを生成し、これらを現世代の集団とする。
- 集団に対し、インスタンスを用いた際の求解時間を評価関数として交叉、選択、突然変異といった処理を行い、次世代の集団を作る。
- ひたすら繰り返せば難しいインスタンスが生成される(はず)
の手法を取ることにした。
遺伝的アルゴリズムの模式図
使用した数理モデル
今回使用する数理モデルとして、多制約ナップザック問題を用いることにした。
定式化は以下の通りとなる。なお、今回の定式化には数理モデリングツールであるJijModelingを使用している。
ソースコード
import jijmodeling as jm
# Problem parameters
p = jm.Placeholder("p", ndim=1, description="Profit values for each item")
w = jm.Placeholder("w", ndim=2, description="Weight matrix where w[m,i] is weight of item i in constraint m")
c = jm.Placeholder("c", ndim=1, description="Capacity constraints for each dimension")
# Get dimensions from placeholders with latex representation
N = p.len_at(0, latex="N") # Number of items from length of profit array
M = w.len_at(0, latex="M") # Number of constraints from first dimension of weight matrix
# Indices
i = jm.Element("i", (0, N), description="Item index")
m = jm.Element("m", (0, M), description="Constraint index")
# Decision variables
x = jm.BinaryVar("x", shape=(N,), description="Binary decision variable for item selection")
# Problem definition
problem = jm.Problem("MultiConstrainedKnapsack", sense=jm.ProblemSense.MAXIMIZE)
# Objective: Maximize total profit
problem += jm.sum(i, p[i] * x[i])
# Constraints: Weight constraints for each dimension
problem += jm.Constraint(
"capacity_constraints",
jm.sum(i, w[m,i] * x[i]) <= c[m],
forall=m
)
problem
ここで、
ランダムにインスタンスを生成した場合
まずは
100変数、10制約の問題を扱ったので、実用上からするとかなり小さい規模の部類のインスタンスであると言える。
* インスタンス
- 変数の数 (N): 100
- 制約の数 (M): 10
* サンプリング
- p_i: 1 ~ 999 からの一様乱数
- w_{m,i}: 1 ~ 999 からの一様乱数
- c_m: 10000 ~ 20000 からの一様変数
- サンプリング回数: 100
* 使用ソルバー
- SCIP (PySCIPOptより使用)
実験コード
import jijmodeling as jm
from ommx_pyscipopt_adapter import instance_to_model, model_to_solution
from ommx.v1 import Instance, Solution
import time
import numpy as np
class Solver:
def __init__(self):
# Problem parameters
p = jm.Placeholder("p", ndim=1, description="Profit values for each item")
w = jm.Placeholder("w", ndim=2, description="Weight matrix where w[m,i] is weight of item i in constraint m")
c = jm.Placeholder("c", ndim=1, description="Capacity constraints for each dimension")
# Get dimensions from placeholders with latex representation
N = p.len_at(0, latex="N") # Number of items from length of profit array
M = w.len_at(0, latex="M") # Number of constraints from first dimension of weight matrix
# Indices
i = jm.Element("i", (0, N), description="Item index")
m = jm.Element("m", (0, M), description="Constraint index")
# Decision variables
x = jm.BinaryVar("x", shape=(N,), description="Binary decision variable for item selection")
# Problem definition
problem = jm.Problem("MultiConstrainedKnapsack", sense=jm.ProblemSense.MAXIMIZE)
# Objective: Maximize total profit
problem += jm.sum(i, p[i] * x[i])
# Constraints: Weight constraints for each dimension
problem += jm.Constraint(
"capacity_constraints",
jm.sum(i, w[m,i] * x[i]) <= c[m],
forall=m
)
self.problem = problem
def solve(self, p, w, c):
# Create instance data
instance_data = {
"p": p,
"w": w,
"c": c
}
# Create interpreter and evaluate problem
interpreter = jm.Interpreter(instance_data)
instance: Instance = interpreter.eval_problem(self.problem)
# PySCIPOpt model に変換
model = instance_to_model(instance)
# PySCIPOptのオプションを指定
model.hideOutput(True)
model.setParam('limits/time', 60)
# 求解
start = time.time()
model.optimize()
elapsed_time = time.time() - start
# 解を取得し、`ommx.v1.Solution`に変換
solution: Solution = model_to_solution(model, instance)
# 経過時間と目的関数値を返す
return elapsed_time, solution.objective
import numpy as np
N = 100
M = 10
p_lower = 1
p_upper = 999
w_lower = 1
w_upper = 999
c_lower = 10000
c_upper = 20000
s = Solver()
elapsed_times = []
for _ in range(100):
instance_data = {
"p": np.random.randint(p_lower, p_upper, size=(N), dtype=np.int32).tolist(),
"w": np.random.randint(w_lower, w_upper, size=(M, N), dtype=np.int32).tolist(),
"c": np.random.randint(c_lower, c_upper, size=(M), dtype=np.int32).tolist(),
}
elapsed_time, _ = s.solve(
instance_data["p"], instance_data["w"], instance_data["c"]
)
elapsed_times.append(elapsed_time)
# Draw Histogram
import matplotlib.pyplot as plt
plt.hist(elapsed_times, bins=50)
plt.xlabel("Elapsed Time [s]")
plt.ylabel("Frequency")
plt.show()
結果は次のようになり、動作環境にはよるものの、筆者の環境では数理最適化ソルバーSCIPを用いておおよそ1秒以内には求まるような問題であることがわかる。
このため、ランダムにインスタンスを生成するのでは比較的容易な問題しか生成されないことがわかる。
ランダムインスタンスの生成結果
遺伝的アルゴリズム
今回は遺伝的アルゴリズムの実装にDEAPと呼ばれるライブラリを用いた。このライブラリを使うと交叉、選択、突然変異といった遺伝的アルゴリズムに必要な操作、およびアップデート処理を自動で行うことができる。
以下の実験環境で動作を行なった。
* 遺伝的アルゴリズム
- 集団数: 100
- 世代数: 40
- 遺伝的操作
- 2点交叉
- 突然変異 (ガウス分布)
- トーナメント選択
- 評価関数
- インスタンスのソルバー (SCIP)を用いた求解時間
- 実行不可能解が出た場合、求解時間として-1を返すことにする
* インスタンス
- 変数の数 (N): 100
- 制約の数 (M): 10
* 初期集団インスタンス
- p_i: 1 ~ 999 からの一様乱数
- w_{m,i}: 1 ~ 999 からの一様乱数
- c_m: 10000 ~ 20000 からの一様変数
実験コード
from deap import base, creator, tools, algorithms
import random
# 行列のサイズ
N = 100
M = 10
matrix_sizes = {
"p": (N, 1), # Nx1 行列
"w": (M, N), # MxN 行列
"c": (M, 1), # Mx1 行列
}
# 個体の1次元ベクトル化のための長さ
p_size = np.prod(matrix_sizes["p"])
w_size = np.prod(matrix_sizes["w"])
c_size = np.prod(matrix_sizes["c"])
total_length = p_size + w_size + c_size
# 評価関数
def evaluate(individual):
# 1次元ベクトルを行列に復元し、インスタンスを生成
p_size = np.prod(matrix_sizes["p"])
w_size = np.prod(matrix_sizes["w"])
c_size = np.prod(matrix_sizes["c"])
p = np.array(individual[:p_size]).reshape(matrix_sizes["p"]).reshape(N).tolist()
w = (
np.array(individual[p_size : p_size + w_size])
.reshape(matrix_sizes["w"])
.tolist()
)
c = (
np.array(individual[p_size + w_size :])
.reshape(matrix_sizes["c"])
.reshape(M)
.tolist()
)
# 問題をソルバーで解き、計測時間を出力する
elapsed_time, _ = s.solve(p, w, c)
return (elapsed_time,)
# DEAPのセットアップ
creator.create("FitnessMax", base.Fitness, weights=(1.0,)) # 難解度を最大化
creator.create("Individual", list, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
# 個体の生成
# 個体生成時の属性をカスタマイズ
def custom_individual():
# p と w を範囲 [1, 999] で生成
p_values = [random.uniform(1, 999) for _ in range(p_size)]
w_values = [random.uniform(1, 999) for _ in range(w_size)]
# c を範囲 [20000, 40000] で生成
c_values = [random.uniform(10000, 20000) for _ in range(c_size)]
return p_values + w_values + c_values
toolbox.register("individual", tools.initIterate, creator.Individual, custom_individual)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
# 遺伝的操作
toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxTwoPoint) # 2点交叉
toolbox.register(
"mutate", tools.mutGaussian, mu=0, sigma=50, indpb=20
) # ガウス分布による突然変異
toolbox.register("select", tools.selTournament, tournsize=3) # トーナメント選択
random.seed(42)
# 初期集団の生成
pop = toolbox.population(n=100)
# 統計データの収集
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", np.mean)
stats.register("std", np.std)
stats.register("min", np.min)
stats.register("max", np.max)
# 最適化の実行
pop, log = algorithms.eaSimple(
pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=40, stats=stats, verbose=True
)
# 最適解の取得
best_ind = tools.selBest(pop, k=1)[0]
print("Best Individual (flattened):", best_ind)
# 最適解を行列に復元
p_size = np.prod(matrix_sizes["p"])
w_size = np.prod(matrix_sizes["w"])
c_size = np.prod(matrix_sizes["c"])
best_p = np.array(best_ind[:p_size]).reshape(matrix_sizes["p"])
best_w = np.array(best_ind[p_size : p_size + w_size]).reshape(matrix_sizes["w"])
best_c = np.array(best_ind[p_size + w_size :]).reshape(matrix_sizes["c"])
print("Best P Matrix:\n", best_p)
print("Best W Matrix:\n", best_w)
print("Best C Matrix:\n", best_c)
print("Best Fitness (difficulty):", best_ind.fitness.values[0])
実験結果 (各世代ごとのインスタンスの、ソルバー求解時間の平均、標準誤差、および最大/最小値)は以下のようになり、確かに世代が増えるごとに求解時間の集団における平均、および最大値も増加していることがわかった。
gen nevals avg std min max
0 100 0.346575 0.3223 0.00823402 1.62108
1 68 0.570045 0.459292 0.066411 1.92539
2 63 0.817366 0.580123 0.0447381 2.78487
3 68 0.98454 0.589473 0.0719602 2.78487
4 65 1.2272 0.780501 0.10636 2.78487
5 55 1.44729 0.905622 0.0994308 3.78067
6 65 1.59616 0.933178 0.0642271 3.78067
7 55 1.88025 1.02517 0.102963 4.81229
8 68 1.70731 1.11041 0.087358 6.30737
9 66 1.94752 1.25047 0.159117 6.30737
10 53 2.32176 1.22699 0.154733 6.30737
...
横軸を世代数、縦軸を求解時間 (秒)としたグラフにすると以下のようになり、明らかに平均、最大値とともにランダムにインスタンスを生成した場合に比べ計算時間が増大しているため、遺伝的アルゴリズムにより求解に時間のかかるインスタンスが生成できていることがわかった。
横軸を世代数としたときの求解時間 (秒)のグラフ
また、筆者の環境でSCIPで20秒ほどかかるインスタンスも生成された。このインスタンスはOMMXを用いて共有しているので、ダウンロードしてPySCIPOptオブジェクトに変換することにより、ローカル環境のSCIPで動作を確認することもできる。
from ommx.artifact import Artifact
from ommx_pyscipopt_adapter import instance_to_model, model_to_solution
# インスタンスをロード
artifact = Artifact.load("ghcr.io/j-i-k-o/2024_adventcalendar/multiconstraint_knapsack_20s:20241221_0928")
# get instance
desc = artifact.get_layer_descriptor(artifact.layers[0].digest)
instance = artifact.get_instance(desc)
# PySCIPOpt model に変換
model = instance_to_model(instance)
# PySCIPOptのオプションを指定
model.hideOutput(False)
model.redirectOutput()
model.setParam("limits/time", 60)
# 求解
start = time.time()
model.optimize()
elapsed_time = time.time() - start
# 解を取得し、`ommx.v1.Solution`に変換
solution = model_to_solution(model, instance)
# Python-MIP modelにも変換が可能であるため、Python-MIPでも求解できる
#from ommx_python_mip_adapter import instance_to_model, model_to_solution
#model = instance_to_model(instance)
#start = time.time()
#model.optimize()
#elapsed_time = time.time() - start
#solution = model_to_solution(model, instance)
まとめ
今回の記事では遺伝的アルゴリズムを用いて恣意的に難しいインスタンスを生成することを試みた。結果として、100変数10制約程度の小さい問題にも関わらずSCIPで20秒ほどかかる問題を簡単に生成することができることを確認できたので、実用上もある程度有用な方法なのではないかと考えている。
ただ、SCIPで解けるからといって任意のソルバーで難しい問題であるとは限らず、例えば上記のOMMXで共有しているインスタンスはCBCを使うと4秒程度で解ける問題である。[2]
とはいえこれに関しても、例えばSCIPの求解時間 + CBCの求解時間のような2つのソルバーを用いた評価指標を用いることである程度解決できるのではと考えている。
また、今回は評価関数として求解時間を選択したが、例えば制約の充足度のような指標を評価関数に用いることで、実行可能解が見つかりにくいインスタンスを恣意的に生成できるかもしれない。
今後やってみると面白そうな内容として、一方は難しいインスタンスを生成するための学習を行い、もう一方はそのインスタンスをより高速に解くための分枝選択を学習するといった敵対的生成ネットワークのような手法を用いることができるのではと考えていたりもする。分枝選択を機械学習で行うことにより性能を向上する研究として V. Nair et al. (2020) のような手法もあるため、ソルバーの性能を向上させる余地は十分にあるように思える。
Discussion