🙄

遺伝的アルゴリズムで難しいインスタンスを生成してみた

2024/12/21に公開

この記事は数理最適化 Advent Calendar 2024の21日目の記事です。

導入

数理最適化の研究や実装において、作成したアルゴリズムの性能検証や異なる手法間のベンチマークを取得するのに使うためのテスト用のインスタンスが欲しい状況が往々にして生じる。

しかしながら、乱数を用いてインスタンスを生成するとしばしば簡単すぎる問題が生成されがちであるため、「そこそこ難しい」 インスタンスを用意するのは意外と厄介なタスクである。

ベンチマークを取得する対象に任意性があるのであればMIPLIB等のインスタンスを用いてベンチマークを取れば良いものの、自作した数理モデルでベンチマークを取る場合には上記の困難は避けられない。

本記事では、遺伝的アルゴリズムを用いて恣意的に難しいインスタンスを生成してみようと思う。

また、今回のコード全てはNotebookとしてGitHub上に公開しており、BinderやColab環境で試すことが可能である。

https://github.com/j-i-k-o/2024_AdventCalendar

インスタンス生成戦略

インスタンスの「難しさ」の定義にはいろいろあるが、一番わかりやすい定義の一つとして「ソルバーの求解時間」が挙げられる。一般に計算時間が長ければ長いほどソルバーが上下界のギャップを埋めるのに苦労しているといえるので、最適解を求めるのが「難しい」インスタンスと言えるだろう。[1]

そのため、

  1. まず乱数を用いてランダムにインスタンスを生成し、これらを現世代の集団とする。
  2. 集団に対し、インスタンスを用いた際の求解時間を評価関数として交叉、選択、突然変異といった処理を行い、次世代の集団を作る。
  3. ひたすら繰り返せば難しいインスタンスが生成される(はず)

の手法を取ることにした。

遺伝的アルゴリズム
遺伝的アルゴリズムの模式図

使用した数理モデル

今回使用する数理モデルとして、多制約ナップザック問題を用いることにした。

定式化は以下の通りとなる。なお、今回の定式化には数理モデリングツールである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
\begin{array}{cccc}\text{Problem:} & \text{MultiConstrainedKnapsack} & & \\& & \max \quad \displaystyle \sum_{i = 0}^{N - 1} p_{i} \cdot x_{i} & \\\text{{s.t.}} & & & \\ & \text{capacity\_constraints} & \displaystyle \sum_{i = 0}^{N - 1} w_{m, i} \cdot x_{i} \leq c_{m} & \forall m \in \left\{0,\ldots,M - 1\right\} \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}& \text{Binary decision variable for item selection}\\\end{array}

ここで、 x_iが0, 1を取る決定変数であり、 p_i, w_{m,i}, c_m が今回決定するべきインスタンスに相当する。それぞれ多制約ナップザック問題の各アイテムの価値、各アイテムの重量、上限値に対応する。

ランダムにインスタンスを生成した場合

まずは p_i, w_{m,i}, c_m をランダムに生成し、どの程度求解に時間を要するかを確認する。以下のような条件で生成を行い、ヒストグラムを表示した。

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) のような手法もあるため、ソルバーの性能を向上させる余地は十分にあるように思える。

脚注
  1. 実用上ではこのようなケースは適当なMIPgapを設定するので、本当であれば「途中で処理を打ち切っても実行可能解が存在していない」といった条件も必要だと思う。実際、今回用いたインスタンスでは実行可能解自体はすぐ見つかるため、MIPgapをソルバーオプションとして設定すれば解決してしまう問題ではある。 ↩︎

  2. 一般にCBCよりもSCIPのほうがパフォーマンスが優れる印象があるので、この結果は個人的に意外だった。 ↩︎

Discussion