Автоматизация поиска гипер-параметров для алгоритма муравьиной колонии

Где-то год назад я столкнулся с задачей по созданию приложения, которое могло бы помочь компании, продающей сладости, строить самый короткий маршрут для логиста, который развозит эти сладости по точкам. Мне нужно было получить координаты точек согласно маршрутного листа и выдать результат как лучше проехать. Где-то месяц я занимался исследованиями. На тот момент я понял немного, фактически только то, что если в маршрутном листе будет больше 10 точек, то с каждой следующей точкой решение задачи в лоб становится нереальным за разумное количество времени. Не только из-за NP-уровня сложности данного типа комбинаторных задач, но и в том числе из-за того, что для получения матрицы расстояний между точками на карте есть всего несколько сервисов с открытым API, но они тоже отдают весьма небольшое количество результатов. После общения с товарищем было придумано решение дробить точки на кластеры и потом искать сначала оптимальный маршрут внутри кластера, а потом уже связывать крайние точки между кластерами. Но пока я занимался исследованиями, задача перестала быть актуальной и эти знания были отложены в сторону.

В следующий раз к задачам на поиск кратчайшего пути я вернулся спустя год. По воле случая пришлось углубиться в современные методы решения подобных задач и я обнаружил целый удивительный мир эвристических, не точных, методов решения логистических задач, о которых год назад и не подозревал. Оказалось, что если решать задачу не точно, а приближенно, то можно получить достаточно близкое к идеальному решение за достаточно разумное количество времени. И в первую очередь меня заинтересовали природные алгоритмы, так как в целом это было достаточно интересно. Как в свое время, когда ученым пришла идея использовать идею передачи сигнала между нейронами головного мозга для создания искусственного интеллекта (хотя, конечно, потом со временем биологи поняли, что мозг работает не так; но математиков это не сильно смутило, так как сама идея нейронные сети искусственного интеллекта (упрощенно будем называть это так, не вдаваясь в детали) прекрасно решали практические задачи), так и здесь, вроде бы есть ощущение, что в природе все сложнее, чем мы можем объяснить простой формулой, однако, это работает. Из природных алгоритмов первым был "алгоритм муравьиной колонии", все остальные крутятся вокруг базовой идеи, заложенной в этом алгоритме. Его и решил исследовать.

Разбирая алгоритм муравьиной колонии дальше, меня заинтересовало то, что этого метод опирается на некоторые параметры, которые можно настраивать. Например, это количество муравьев, которые принимают участие в поиске маршрута, а также количество итераций. И уже существует отдельное научное направление по поиску методов, с помощью которых можно находить оптимальные параметры для данного алгоритма.

Не буду углубляться в объяснение самого алгоритма, об этом вы можете прочитать в Википедии.

Мне пришла мысль использовать методы для поиска оптимальных гипер-параметров из классического машинного обучения в связке с алгоритмом муравьиной колонии. Эта идея в легла в основу моей диссертации. И сейчас я хочу поделиться с вами своими наблюдениями.

Для начала я решил остановиться на давно известных методах - Grid Search, Random Search и Bayesian Optimization. Опять же, кто хочет изучить детальнее, как работает каждый из этих методов, в интернете много на этот счет информации. Не хочу перегружать статью. Например, можете прочитать здесь. Ранее исследователями были проведены замеры скорости и точности работы этих трех алгоритмов, в результате которых более предпочительным для использования оказался Bayesian Optimization. Его я и решил применить.

Чтобы проверить свою идею как можно быстрее, я решил воспользоваться уже готовой библиотекой hyperopt, которая позволяет использовать данный метод для поиска гипер-параметров.

В итоге у меня родился такой код:

NUM_POINTS = 500
n = NUM_POINTS

# Генерируем тестовые данные
import numpy as np
from scipy import spatial


def save_to_file(list, file_name):
    # print(list, type(list))
    np.savetxt(f'{ROOT_DIR}/{file_name}_{NUM_POINTS}.csv', list, delimiter=",")


def get_points_coordinate(num_points=NUM_POINTS):
    # num_points - количество вершин
    points_coordinate = np.random.rand(num_points, 2)  # генерация рандомных вершин

    return points_coordinate


def get_distance_matrix(points_coordinate):

    # вычисление матрицы расстояний между вершин
    distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')

    return distance_matrix


POINTS_COORDINATE = get_points_coordinate(NUM_POINTS)

save_to_file(POINTS_COORDINATE, 'POINTS_COORDINATE')

DISTANCE_MATRIX = get_distance_matrix(POINTS_COORDINATE)

save_to_file(DISTANCE_MATRIX, 'DISTANCE_MATRIX')

# Реализация алгоритма
import time
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from hyperopt import hp, fmin, tpe, Trials, STATUS_OK
from numpy import genfromtxt


class ACO_TSP:  # класс алгоритма муравьиной колонии для решения задачи коммивояжёра
    """
    Код написан в соответствие с описанием алгоритма и формулами выше. В алгоритме есть два важных фактора:
    феромон τ и видимость η. Феромон τ относится к оставшейся информации на каждом пути, по которому проходят муравьи.
    Видимость η обозначает обратное расстояние между узлами. Вероятность P состоит из произведения
    феромонаτ и видимости η и их суммы. При кратчайшем пути, высоком уровне феромона τ, а также высокой
    видимости η вероятность P становится выше. Следовательно, в строках 39-40 происходит подсчёт вероятности P.

    Затем, строка 41 выбирает последовательность кратчайших путей при помощи вероятности P.

    Наконец, строки 55-64 обновляют феромон τ с коэффициентом испарения ρ.

    Параметр Q – это сумма феромонов в постоянном значении. Вернувшись к строке 37,
    новая последовательность маршрутов определена для следующего вычисления.

    """

    def __init__(self, func, n_dim, size_pop=10, max_iter=20, distance_matrix=None, alpha=1, beta=2, rho=0.1):
        self.func = func
        self.n_dim = n_dim  # количество городов
        self.size_pop = size_pop  # количество муравьёв
        self.max_iter = max_iter  # количество итераций
        self.alpha = alpha  # коэффициент важности феромонов в выборе пути
        self.beta = beta  # коэффициент значимости расстояния
        self.rho = rho  # скорость испарения феромонов
        self.distance_matrix = distance_matrix  # матрица расстояний

        self.prob_matrix_distance = 1 / (distance_matrix + 1e-10 * np.eye(n_dim, n_dim))

        # Матрица феромонов, обновляющаяся каждую итерацию
        self.Tau = np.ones((n_dim, n_dim))
        # Путь каждого муравья в определённом поколении
        self.Table = np.zeros((size_pop, n_dim)).astype(int)
        self.y = None  # Общее расстояние пути муравья в определённом поколении
        self.generation_best_X, self.generation_best_Y = [], []  # фиксирование лучших поколений
        self.x_best_history, self.y_best_history = self.generation_best_X, self.generation_best_Y
        self.best_x, self.best_y = None, None
        self.iter_distance = []  # итерационное расстояние муравья в определённом поколении
        self.distance_dynamic = []
        self.score = []

    def run(self, max_iter=None):
        self.max_iter = max_iter or self.max_iter
        for i in range(self.max_iter):
            # вероятность перехода без нормализации
            prob_matrix = (self.Tau ** self.alpha) * (self.prob_matrix_distance) ** self.beta
            for j in range(self.size_pop):  # для каждого муравья
                # точка начала пути (она может быть случайной, это не имеет значения)
                self.Table[j, 0] = 0
                for k in range(self.n_dim - 1):  # каждая вершина, которую проходят муравьи
                    # точка, которая была пройдена и не может быть пройдена повторно
                    taboo_set = set(self.Table[j, :k + 1])
                    # список разрешённых вершин, из которых будет происходить выбор
                    allow_list = list(set(range(self.n_dim)) - taboo_set)
                    prob = prob_matrix[self.Table[j, k], allow_list]
                    prob = prob / prob.sum()  # нормализация вероятности
                    next_point = np.random.choice(allow_list, size=1, p=prob)[0]
                    self.Table[j, k + 1] = next_point

            # рассчёт расстояния
            y = np.array([self.func(i, self.distance_matrix) for i in self.Table])

            # фиксация лучшего решения
            index_best = y.argmin()
            x_best, y_best = self.Table[index_best, :].copy(), y[index_best].copy()

            # добавляем найденную лучшую дистанцию на этой итерации в список для графика
            self.iter_distance.append(y_best)
            # производим оценку качества найденного оптимума (лучше, чем придыдущий или нет)
            if y_best < get_min_int_in_list(self.distance_dynamic, y_best):
                self.distance_dynamic.append(y_best)
            else:
                self.distance_dynamic.append(self.distance_dynamic[-1])
            if len(self.distance_dynamic) > 1:
                self.score.append(self.distance_dynamic[-1] - self.distance_dynamic[-2])
            else:
                self.score.append(1)

            self.generation_best_X.append(x_best)
            self.generation_best_Y.append(y_best)

            # подсчёт феромона, который будет добавлен к ребру
            delta_tau = np.zeros((self.n_dim, self.n_dim))
            for j in range(self.size_pop):  # для каждого муравья
                for k in range(self.n_dim - 1):  # для каждой вершины
                    # муравьи перебираются из вершины n1 в вершину n2
                    n1, n2 = self.Table[j, k], self.Table[j, k + 1]
                    delta_tau[n1, n2] += 1 / y[j]  # нанесение феромона
                # муравьи ползут от последней вершины обратно к первой
                n1, n2 = self.Table[j, self.n_dim - 1], self.Table[j, 0]
                delta_tau[n1, n2] += 1 / y[j]  # нанесение феромона

            self.Tau = (1 - self.rho) * self.Tau + delta_tau

        best_generation = np.array(self.generation_best_Y).argmin()
        self.best_x = self.generation_best_X[best_generation]
        self.best_y = self.generation_best_Y[best_generation]

        score = (self.score.count(0) - len(self.score)) / len(self.score)

        total_distance = self.get_total_distance()

        return self.best_x, self.best_y, self.iter_distance, self.y_best_history, score, total_distance

    def get_total_distance(self):
        # А теперь посчитаем саму длину маршрута
        best_points_ = np.concatenate([self.best_x, [self.best_x[0]]])
        best_points_coordinate = POINTS_COORDINATE[best_points_, :]

        # Построим оптимальный маршрут в виде списка элементов в массиве точек для дальнейшнего просчете расстояния оптимального маршрута
        list_of_points_coordinate = list(POINTS_COORDINATE.tolist())
        list_of_best_points_coordinate = list(best_points_coordinate.tolist())

        list_of_index_set_of_points = []
        for el in list_of_best_points_coordinate:
            list_of_index_set_of_points.append(list_of_points_coordinate.index(el))

        total_distance = 0
        for i in range(1, len(list_of_index_set_of_points)):
            total_distance += cal_euclidean_distance_between_two_points(best_points_coordinate[i],
                                                                        best_points_coordinate[i - 1])

        return total_distance

    fit = run


# вычисление длины пути
def cal_total_distance(routine, distance_matrix):
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])


def cal_euclidean_distance_between_two_points(point_1, point_2):
    return np.sqrt(sum(pow(a - b, 2) for a, b in zip(point_1, point_2)))


def get_min_int_in_list(list, y_best):
    if len(list) < 1:
        return y_best + 1
    else:
        return min(list)


def aco_tsp_basic(num_points, distance_matrix, size_pop=3, max_iter=50):
    # создание объекта алгоритма муравьиной колонии
    aca = ACO_TSP(func=cal_total_distance,
                  n_dim=num_points,
                  size_pop=size_pop,  # количество муравьёв
                  max_iter=max_iter,
                  distance_matrix=distance_matrix)
    best_x, best_y, iter_distance, y_best_history, score, total_distance = aca.run()

    return aca, best_x, best_y, iter_distance, y_best_history, score, total_distance


def aco_tsp_opt(params, num_points=NUM_POINTS, distance_matrix=DISTANCE_MATRIX):
    print(params)
    # создание объекта алгоритма муравьиной колонии
    aca = ACO_TSP(func=cal_total_distance,
                  n_dim=num_points,
                  distance_matrix=distance_matrix,
                  **params)
    best_x, best_y, iter_distance, y_best_history, score, total_distance = aca.run()

    return {'loss': total_distance, 'params': params, 'status': STATUS_OK}


def plot(aca, best_x, iter_distance, y_best_history):
    # Вывод результатов на экран
    fig, ax = plt.subplots(1, 3)
    best_points_ = np.concatenate([best_x, [best_x[0]]])
    best_points_coordinate = POINTS_COORDINATE[best_points_, :]

    # Выведем на графике данные о точках и как они соединялись между собой с точки зрения оптимального прохождения маршрута
    for index in range(0, len(best_points_)):
        ax[0].annotate(best_points_[index], (best_points_coordinate[index, 0], best_points_coordinate[index, 1]))
    ax[0].plot(best_points_coordinate[:, 0],
               best_points_coordinate[:, 1], 'o-r')

    # Выведем куммулятивный минимум расстояний, которые находили за все время работы
    pd.DataFrame(aca.y_best_history).cummin().plot(ax=ax[1], title="Кумулятивно мінімально/оптимальна відстань",
                                                   legend=False)

    # Выведем статистику расстояний, которые находили за все время работы
    pd.DataFrame(iter_distance).plot(ax=ax[2], title="Статистика відстаней", legend=False)

    # изменение размера графиков
    plt.rcParams['figure.figsize'] = [50, 50]
    plt.show()


start_time = time.time()  # сохранение времени начала выполнения

search_space = {
    'size_pop': hp.choice(label='size_pop',
                          options=np.arange(1, 251, 5, dtype=int)),
    'max_iter': hp.choice(label='max_iter',
                          options=np.arange(1, 251, 5, dtype=int)),
}


# запускаем hyperopt
trials = Trials()
best = fmin(
    # функция для оптимизации
    fn=aco_tsp_opt,
    # пространство поиска гиперпараметров
    space=search_space,
    # алгоритм поиска
    algo=tpe.suggest,
    # число итераций
    # (можно ещё указать и время поиска)
    max_evals=20,
    # куда сохранять историю поиска
    trials=trials,
    # progressbar
    show_progressbar=True
)


def df_results(hp_results):
    """
    Отображаем результаты hyperopt в формате DataFrame

    :hp_results: результаты hyperop
    :return: pandas DataFrame
    """
    results = pd.DataFrame([{**x, **x['params']} for x in  hp_results])
    results.sort_values(by=['loss'], ascending=True, inplace=True)
    return results


results = df_results(trials.results)
print(results)
print("time of execution: %s seconds" % abs(time.time() - start_time))  # вычисление времени выполнения

Сам код для алгоритма муравьиной колонии мне было лень писать, поэтому решил воспользоваться уже готовым решением. Мне было интересно просто проверить свою гипотезу.

Чтобы понять насколько данное направление имеет в целом смысл, решил сравнить полученные результаты с тем, что мне выдал бы "жадный" алгоритм - метод ближайшего соседа. Он считается достаточно точным методом, еще и работающим достаточно быстро относительно других.

Получился такой результат:

Данные просчеты проводились с помощью бесплатной версии Google Colab.

Что мы видим в этой таблице?

Во-первых, смысл в подборе параметров для АМК есть. Неоптимизированный вариант дает очень плохой результат, а оптимизированный - даже лучше, чем метод ближайшего соседа.

Во-вторых, подбирать параметры самому тяжело - действительно Bayesian Optimization помогает автоматизировать процесс подбора оптимальных гипер-параметров.

В-третьих, есть и ряд проблем - оптимизированный АМК (как правило, это большее кол-во итераций или особей) работает намного дольше, чем метод ближайшего соседа. Но уже хорошо известно, что АМК может работать параллельно (то есть особи муравьев могут запускаться параллельно, не последовательно) и это значительно ускоряет получение результата. Для параллельных вычислений нам могут понадобиться методы оптимизации (поиска оптимальных параметров) более продвинутые, чем Bayesian Optimization (потому что этот метод исходя из своей логики опирается на предыдущие поиски), хотя, мы можем применить Bayesian Optimization и для параллельных исчислений тоже (в этом смысле у меня появились определенные идеи).

Поэтому на этом исследование не прекращается - еще есть куда развивать идею дальше.

Буду делиться в последующем полученными результатами.