メタヒューリスティクス開発者向けのツールopytimizerを使ってみる

本記事について

メタヒューリスティクスを開発する際に、開発コストを下げるツールがないかと探していたら、pythonパッケージであるopytimizerを見つけました。簡単に触ってみて面白いなと思ったので紹介します。０から全て開発するよりもopytimizerを使った方が少し便利かなという印象です。
本記事は、pythonに慣れていて、メタヒューリスティクス等の数理最適化アルゴリズムの開発者やアルゴリズム実装に関心がある方向けに書いています。

はじめに

メタヒューリスティクスとは、数理最適化用語集によると以下のように定義されています。具体的にはシミュレーテッドアニーリングやタブーサーチ等が挙げられます。

ヒューリスティクスとは実験的手法によって解を導出する方法であり，普通ある問題に特化した方法をいう．それに対しメタヒューリスティクスは同様に実験的手法ではあるが，その枠組みを拡張して様々な問題に応用できる汎用的な近似アルゴリズムのことをいう．

実務の課題を数理最適化の技術を使って解決しようとすると、既存の手法では速度が遅い、実行可能解すら見つからない...といったケースに直面するときがあります。性能を受け入れて運用する時もあれば、ヒューリスティクスやメタヒューリスティクスを実装して性能を改善するときもあります。

メタヒューリスティクスの開発では開発コストがかかる点がボトルネックになってきます。手法の良し悪しは、実際に最適化問題を解いてみないとわからないこともよくあり、解けない原因を調査して改善...を繰り返していくので、難しい問題だと開発量が膨れ上がります。限られた期間内で性能を出すためには、もちろん開発スキルを向上させることも重要ですが、開発コストを削減できるかどうかも重要になってきます。

opytimizerとは

コンセプト

開発コストを削減できそうなツールとして、pythonパッケージのopytimizer（読み方はわからない...）というものがあります。README.mdを読むとメタヒューリスティクス開発者向けのツールであることが言及されています。

Did you ever reach a bottleneck in your computational experiments? Are you tired of selecting suitable parameters for a chosen technique? If yes, Opytimizer is the real deal! This package provides an easy-to-go implementation of meta-heuristic optimizations. From agents to search space, from internal functions to external communication, we will foster all research related to optimizing stuff.
Use Opytimizer if you need a library or wish to:
Create your optimization algorithm;
Design or use pre-loaded optimization tasks;
Mix-and-match different strategies to solve your problem;
Because it is fun to optimize things.

用語

実装する上で見慣れない用語が出てきたので私の解釈を記載しておきます。

• Agent: おそらく解の情報。n_agentsは解プールに保存する解の個数
• fit: 目的関数値
• position: 解（$x=1$などの情報）

このツールを使ってできること

これだけだと具体的に何ができるのかイメージがわかなかったので、opytimizerを使って、TSP (巡回セールスマン問題)を例に2-opt近傍を実装してみました。特徴としては以下の通りです。

• 目的関数を自由に定義可能（解を入力として計算式を実装するだけで良い）
• 制約条件をboolianで判定する形で定義可能（制約を満たすならTrue、そうでないならFalse）
• 探索空間を定義可能
• 解を更新する処理を実装可能

上記をユーザーで定義したら、あとはopytimizerに用意されているstart関数を呼べば最適化が実行される。

使い方

TSPに対する2-opt近傍の実装を例に説明します。2-optのイメージはこちらの記事等を参照してください。

インストール方法

pip install opytimizer

実装例

データセット

今回は こちらのコードを利用させていただきます。

from geopy.geocoders import Nominatim
import numpy as np

geolocator = Nominatim(user_agent="python")

# set the name list of traveling points
points = ['茨城県', '栃木県', '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県']

# get the latitude and longitude
latlng_list = []
for point in points:
    location = geolocator.geocode(point)
    if location is not None:
        latlng_list.append([ location.latitude, location.longitude ])
# make distance matrix
num_points = len(points)
inst_d = np.zeros((num_points, num_points))
for i in range(num_points):
    for j in range(num_points):
        a = np.array(latlng_list[i])
        b = np.array(latlng_list[j])
        inst_d[i][j] = np.linalg.norm(a-b)
# normalize distance matrix
inst_d = (inst_d-inst_d.min()) / (inst_d.max()-inst_d.min())

geo_data = {'points': points, 'latlng_list': latlng_list}
instance_data = {'d': inst_d}

var_map = {(i, t): num_points * i + t for i in range(num_points) for t in range(num_points)}

目的関数の定義

総移動距離の最小化を実装します。基本的には、解x (np.array型)を入力し、目的関数の計算結果 (float型)を返す関数を定義すれば良いです。今回は、距離マトリックスを持たせたかったので、以下のようにクラスを定義して、特殊メソッド__call__を実装することで関数としても呼べるようにしています。注意点として、xのフォーマットに触れておきます。xは、例えば np.array([[0], [1]])のような多次元配列になっています。それぞれの要素を$x^1, x^2$と仮に置くと、１つ目の要素は$x^1$の次元数分の値が入ります。２つ目の要素も同様です。$x^1_0, x^1_1, ...., x^1_n, x^2_0, x^2_1, ...x^2_n$で各要素は次元数$n$、配列の長さは変数の数（この例だと$x^1,x^2$の２つ）というフォーマットのようです。今回は、$x_{it}$を$x^1, x^2, ...$という感じで割り当てました。
※筆者はフォーマットの理解を間違えて実装してしまってハマりました...

class Distance:

    def __init__(self, distance_matrix: np.array, var_map: dict, num_nodes: int):
        self.distance_matrix = distance_matrix
        self.var_map = var_map
        self.num_nodes = num_nodes

    def __call__(self, x: np.array) -> float:
        travel_distance = 0.0
        route = []
        for p, values in enumerate(x):
            if values[0] == 0:
                # n_dimension = 1
                continue
            city = p % self.num_nodes
            route.append(city)
        if len(route) != self.num_nodes:
            return 100000
        for i, city in enumerate(route):
            next_city = route[(i+1)%self.num_nodes]
            travel_distance += self.distance_matrix[city][next_city]
        return travel_distance

制約条件の定義

以下の制約条件を定義しています。

• 各時刻で訪問できる都市は１つだけ
• 各都市に訪問するのは１回だけ

目的関数と同様に、基本的には解x (np.array型)を入力し、制約を満たすならTrue、違反するならFalseを返す関数を実装すれば良いです。今回も追加で情報を持たせたかったので以下のように実装しました。

class TimeDirectionOneHot:

    def __init__(self, var_map: dict, num_nodes: int):
        self.var_map = var_map
        self.num_nodes = num_nodes

    def __call__(self, x: np.array) -> bool:
        for i in range(self.num_nodes):
            if np.sum([x[self.var_map[(i, t)]][0] for t in range(self.num_nodes)]) != 1:
                return False
        return True


class NodeDirectionOneHot:

    def __init__(self, var_map: dict, num_nodes: int):
        self.var_map = var_map
        self.num_nodes = num_nodes

    def __call__(self, x: np.array) -> bool:
        for t in range(self.num_nodes):
            if np.sum([x[self.var_map[(i, t)]][0] for i in range(self.num_nodes)]) != 1:
                return False
        return True

探索空間の定義

探索空間を定義していきます。opytimizerに用意されているクラス（BooleanSpaceなど）を使用するのもできますが、今回は初期解の生成を実装したかったので自分で定義しました。初期解の生成は、Spaceクラスにある_initialize_agentsをオーバーライドする形で実装できます。処理のしやすさで、routesという属性値を持たせています。初期解の生成方法は、訪問する順番にランダムで都市を選んでいくというシンプルな設計にしました。

import copy
from typing import List, Optional
import random

import numpy as np

from opytimizer.core import Space, Function
from opytimizer.utils import logging

logger = logging.get_logger(__name__)


class TourSpace(Space):
    """A TourSpace class for agents, variables and methods
    related to the boolean search space.

    """

    def __init__(
        self, n_agents: int, n_variables: int, n_cities: int,
        objective_function: Function,
        mapping: Optional[List[str]] = None
    ) -> None:
        """Initialization method.

        Args:
            n_agents: Number of agents.
            n_variables: Number of decision variables.
            mapping: String-based identifiers for mapping variables' names.

        """

        logger.info("Overriding class: Space -> TourSpace.")

        n_dimensions = 1
        lower_bound = np.zeros(n_variables)
        upper_bound = np.ones(n_variables)
        self.n_cities = n_cities
        self.objective_function = objective_function
        self.routes = []

        super(TourSpace, self).__init__(
            n_agents, n_variables, n_dimensions, lower_bound, upper_bound, mapping
        )

        self.build()

        logger.info("Class overrided.")

    def _initialize_agents(self) -> None:
        """Initializes agents with their positions and defines a best agent."""

        for agent in self.agents:
            for j in range(self.n_variables):
                agent.position[j] = [0]

        for agent in self.agents:
            added_cities = set()
            route = []
            for n in range(self.n_cities):
                city = random.randint(0, self.n_cities-1)
                while city in added_cities:
                    city = random.randint(0, self.n_cities-1)
                index = n * self.n_cities + city
                agent.position[index][0] = 1
                agent.fit = self.objective_function(agent.position)
                added_cities.add(city)
                route.append(city)
            self.routes.append(route)

        argmin = np.argmin([agent.fit for agent in self.agents])
        self.best_agent = copy.deepcopy(self.agents[argmin])

2-optの実装

Optimizerというクラスを継承し、update関数を実装するだけで基本的には大丈夫です。現在の解の情報はspaceのagentsに入っています。2-opt適用後の解に更新すれば良いです。2-optは、各々の解に対して、ランダムで辺を選択、2-optが適用できる辺を順に調べて一番コストが下がる辺を選択して2-optを適用する、という流れで実装しました。

import copy
from typing import Any, Dict, Optional

import numpy as np

import opytimizer.math.random as r
from opytimizer.core import Optimizer
from opytimizer.core.function import Function
from opytimizer.core.space import Space
from opytimizer.utils import logging


class TwoOpt(Optimizer):
    """A TwoOpt class, inherited from Optimizer.

    This is the designed class to define TwoOpt-related
    variables and methods.

    """

    def __init__(self, params: Optional[Dict[str, Any]] = None) -> None:
        """Initialization method.

        Args:
            params: Contains key-value parameters to the meta-heuristics.

        """

        logger.info("Overriding class: Optimizer -> TwoOpt.")

        super(TwoOpt, self).__init__()

        self.build(params)

        logger.info("Class overrided.")

    def update(self, space: TourSpace, function: Function) -> None:
        """Wraps Passing Vehicle Search over all agents and variables.

        Args:
            space: Space containing agents and update-related information.
            function: A Function object that will be used as the objective function.

        """

        d_mat = function.pointer.distance_matrix
        for i, agent in enumerate(space.agents):
            route = space.routes[i]
            visit_pos1 = random.randint(0, space.n_cities-1)
            if visit_pos1 == space.n_cities-1:
                edge1 = (route[visit_pos1], route[0])
            else:
                edge1 = (route[visit_pos1], route[visit_pos1+1])

            best_reduction = 0.0
            best_route = copy.deepcopy(route)
            for visit_pos2 in range(space.n_cities):
                city = route[visit_pos2]
                if city in edge1:
                    continue
                next_city = route[0] if visit_pos2 == space.n_cities - 1 else route[visit_pos2+1]
                if next_city in edge1:
                    continue
                edge2 = (city, next_city)
                if visit_pos1 < visit_pos2:
                    first, second = (edge1, edge2)
                    first_pos, second_pos = visit_pos1, visit_pos2
                else:
                    first, second = (edge2, edge1)
                    first_pos, second_pos = visit_pos2, visit_pos1

                reduction = (d_mat[first[0]][second[0]] + d_mat[first[1]][second[1]]
                             - d_mat[first[0]][first[1]] - d_mat[second[0]][second[1]])
                logger.info(f"{first}, {second}")
                logger.info(f"Before {d_mat[first[0]][first[1]] + d_mat[second[0]][second[1]]}")
                logger.info(f"After {d_mat[first[0]][second[0]] + d_mat[first[1]][second[1]]}")
                if reduction < best_reduction:
                    best_route = (
                        route[:first_pos+1] + [second[0]] + list(reversed(route[first_pos+2:second_pos]))
                        + [first[1]] + route[second_pos+1:]
                    )
                    best_reduction = reduction

            if best_reduction < 0.0:
                logger.info(f"Update solution for Agent {i}!")
                logger.info(f"route: {route} -> {best_route}")
                space.routes[i] = copy.deepcopy(best_route)
                old_fit = agent.fit
                position = np.zeros((space.n_variables, 1))
                for n, city in enumerate(space.routes[i]):
                    index = space.n_cities * n + city
                    position[index] = [1]
                agent.position = copy.deepcopy(position)
                agent.fit = function(agent.position)
                logger.info(f"fit: {old_fit} -> {agent.fit}")

最適化の実行

以下のような形で実行できます。実行するとログが自動で生成されます。

import numpy as np

from opytimizer import Opytimizer
from opytimizer.functions import ConstrainedFunction
from opytimizer.optimizers.boolean import BMRFO, BPSO, UMDA
from opytimizer.spaces import BooleanSpace

# Objective function
cost = Distance(distance_matrix=inst_d, var_map=var_map, num_nodes=num_points)

# Defines a constraint function that returns a boolean
# whether the constraint is valid or not
onhot1 = TimeDirectionOneHot(var_map=var_map, num_nodes=num_points)
onehot2 = NodeDirectionOneHot(var_map=var_map, num_nodes=num_points)


# Random seed for experimental consistency
np.random.seed(0)

# Number of agents and decision variables
n_agents = 5
n_variables = num_points ** 2

# Creates the space, optimizer and function
space = TourSpace(n_agents, n_variables, num_points, cost)
optimizer = TwoOpt()
function = ConstrainedFunction(cost, [onhot1, onehot2], penalty=100.0)

# Bundles every piece into Opytimizer class
opt = Opytimizer(space, optimizer, function, save_agents=False)

# Runs the optimization task
opt.start(n_iterations=10)

# Best solution
print(opt.space.best_agent.fit)
print(opt.space.best_agent.position)

以下のようなログが表示されます。標準では、各反復の目的関数の値や解が表示されます。追加で、目的関数の変化もログに出力しています。目的関数値が約4.64から約2.87まで落ちました。

2025-01-19 20:15:25,411 - __main__ — INFO — Overriding class: Space -> TourSpace.
2025-01-19 20:15:25,414 - opytimizer.core.space — DEBUG — Agents: 5 | Size: (49, 1) | Lower Bound: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0.] | Upper Bound: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1.] | Mapping: ['x0', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9', 'x10', 'x11', 'x12', 'x13', 'x14', 'x15', 'x16', 'x17', 'x18', 'x19', 'x20', 'x21', 'x22', 'x23', 'x24', 'x25', 'x26', 'x27', 'x28', 'x29', 'x30', 'x31', 'x32', 'x33', 'x34', 'x35', 'x36', 'x37', 'x38', 'x39', 'x40', 'x41', 'x42', 'x43', 'x44', 'x45', 'x46', 'x47', 'x48'] | Built: True.
2025-01-19 20:15:25,416 - __main__ — INFO — Class overrided.
2025-01-19 20:15:25,417 - __main__ — INFO — Overriding class: Optimizer -> TwoOpt.
2025-01-19 20:15:25,417 - opytimizer.core.optimizer — DEBUG — Algorithm: TwoOpt | Custom Parameters: None | Built: True.
2025-01-19 20:15:25,417 - __main__ — INFO — Class overrided.
2025-01-19 20:15:25,418 - opytimizer.functions.constrained — INFO — Overriding class: Function -> ConstrainedFunction.
2025-01-19 20:15:25,418 - opytimizer.core.function — INFO — Creating class: Function.
2025-01-19 20:15:25,419 - opytimizer.core.function — DEBUG — Function: Distance | Built: True.
2025-01-19 20:15:25,419 - opytimizer.core.function — INFO — Class created.
2025-01-19 20:15:25,420 - opytimizer.functions.constrained — DEBUG — Constraints: [<__main__.TimeDirectionOneHot object at 0x119199750>, <__main__.NodeDirectionOneHot object at 0x11919a610>] | Penalty: 100.0.
2025-01-19 20:15:25,420 - opytimizer.functions.constrained — INFO — Class overrided.
2025-01-19 20:15:25,420 - opytimizer.opytimizer — INFO — Creating class: Opytimizer.
2025-01-19 20:15:25,421 - opytimizer.opytimizer — DEBUG — Space: <__main__.TourSpace object at 0x1191b2dd0> | Optimizer: <__main__.TwoOpt object at 0x1191e1f90>| Function: <opytimizer.functions.constrained.ConstrainedFunction object at 0x1191e8a90>.
2025-01-19 20:15:25,421 - opytimizer.opytimizer — INFO — Class created.
2025-01-19 20:15:25,422 - opytimizer.opytimizer — INFO — Starting optimization task.
  0%|          | 0/10 [00:00<?, ?it/s]

(中略)

2025-01-19 20:15:25,456 - __main__ — INFO — fit: 4.642390100592843 -> 4.249882906514916

（中略）

100%|##########| 10/10 [00:00<00:00, 47.03it/s, fitness=2.87]
2025-01-19 20:15:25,640 - opytimizer.opytimizer — INFO — Optimization task ended.
2025-01-19 20:15:25,641 - opytimizer.opytimizer — INFO — It took 0.21631526947021484 seconds.

さいごに

opytimizerを使ってみて面白いなと思いましたし、以下の点は良いと思いました。

• 数式を書かなくても最適化ができる
• 一部だけ実装すれば最適化を実行できる
• ログを自動で出力してくれる

今回は自分で実装してみましたがいくつか実装されているアルゴリズムがあるようなのでそちらを使ってみても良いかもしれません。速度を追求するならCはRust等の言語で実装した方が良い、という話はありますが、pythonに慣れている方は試してみてください。