Практика 13. Эвристические алгоритмы. Метод отжига. Муравьиный алгоритм

1 Эвристические алгоритмы

Эвристический алгоритм — это алгоритм решения задачи, правильность которого для всех возможных случаев не доказана, но про который известно, что он даёт достаточно хорошее решение в большинстве случаев. Эвристические алгоритмы обычно применяются при решении плохо формализованных или сложных задач.

С направленными методами поиска решения курс уже встречался: жадные алгоритмы (практика 12) на каждом шаге выбирают локально лучший вариант, а метод градиентного спуска (практика 4) шагает против градиента целевой функции. Оба подхода разделяют общую слабость — они останавливаются в первом же локально-хорошем решении, которое не обязано быть глобально лучшим. Сегодняшняя практика посвящена двум эвристическим алгоритмам, которые специально устроены так, чтобы из локальных «ям» выбираться: методу имитации отжига и муравьиному алгоритму.

2 Алгоритм имитации отжига

Алгоритм имитации отжига (англ. simulated annealing) — эвристический алгоритм глобальной оптимизации, особенно эффективный при решении дискретных и комбинаторных задач.

Алгоритм вдохновлён процессом отжига в металлургии — техники, заключающейся в нагревании и постепенном охлаждении металла для увеличения его кристаллизованности и уменьшения дефектов.

Симулирование отжига в переборных задачах может быть использовано для приближённого нахождения глобального минимума функций с большим количеством свободных переменных.

Алгоритм вероятностный и не даёт почти никаких гарантий сходимости, однако хорошо работает на практике при решении NP-полных задач. Иногда на контестах им удаётся сдать сложные комбинаторные задачи, у которых есть нормальное решение.

Рассуждения о методе

Как и всё гениальное, данный метод подсмотрен у матушки природы. За основу взят процесс кристаллизации вещества, который, в свою очередь, «приручили» хитрые металлурги для повышения однородности металла.

Напомню, что у каждого металла есть кристаллическая решетка — если совсем коротко, она описывает геометрическое положение атомов вещества. Совокупность позиций всех атомов будем называть состоянием системы, каждому состоянию соответствует определённый уровень энергии. Цель отжига – привести систему в состояние с наименьшей энергией. Чем ниже уровень энергии, тем «лучше» кристаллическая решетка, т.е. тем меньше у нее дефектов и прочнее металл.

В ходе «отжига» металл сначала нагревают до некоторой температуры, что заставляет атомы кристаллической решетки покинуть свои позиции. Затем начинается медленное и контролируемое охлаждение. Атомы стремятся попасть в состояние с меньшей энергией, однако, с определённой вероятностью они могут перейти и в состояние с большей. Эта вероятность уменьшается вместе с температурой. Переход в худшее состояние, как ни странно, помогает в итоге отыскать состояние с энергией меньшей, чем начальная. Процесс завершается, когда температура падает до заранее заданного значения.

Как видите, в рамках данного процесса происходит минимизация энергии. Меняем энергию на нашу целевую функцию, и voilà! Или нет? Давайте задумаемся, а чем же так хорош столь мудрёный процесс? Почему бы нам, например, всегда не переходить в состояния с меньшей энергией?

Примерно так работает имеющий широкое распространение метод градиентного спуска (практика 4).

На просторах интернета мне удалось отыскать совершенно чудесную картинку, которая как нельзя лучше иллюстрирует все плюсы алгоритма.

Давайте на минуточку представим себя лыжником, несущемся по крутому склону нашей целевой функции. Задача спуститься в самую нижнюю точку (глобальный минимум). Но, как видите, склон нам попался не совсем ровный, есть на нём низменности, заехав в которые, нам может показаться, что мы уже у подножья. Дабы не остаться в дураках, нам придётся проявить любопытство и вскарабкаться на небольшой холмик, чтобы продолжить спуск.

Чтобы метафора стала окончательно понятна, в примере с лыжником склон – функция измерения энергии вещества (она же целевая функция). А вероятностный переход в состояние с большей энергией – карабканье на холм справа от низменности.

Как видите, такой алгоритм позволяет нам с большой вероятностью избежать «застревания» в локальном минимуме.

Общее описание алгоритма

Пусть имеется некоторая функция f(x) от состояния x, которую требуется минимизировать. В качестве базового решения выбирается некоторое состояние x₀ — обычно случайное. Вводится температура t — положительное действительное число, которое будет уменьшаться по мере оптимизации и определять вероятность перехода в худшее соседнее состояние.

Алгоритм выполняется до тех пор, пока не будет достигнуто приемлемое решение или не закончится отведённое время. На каждой итерации температура уменьшается по некоторой функции охлаждения, например tₖ = γ ⋅ tₖ₋₁, где γ — число, близкое к единице (скажем, 0.99), или tₖ = t₀ / (1 + k). Затем выбирается случайное соседнее состояние y, получаемое из текущего x небольшим изменением. Если y не хуже (f(y) ≤ f(x)), переход выполняется всегда. В противном случае переход осуществляется лишь с некоторой вероятностью — тем меньшей, чем хуже новое состояние и чем ниже температура.

Ещё раз, суть процесса отжига заключается в следующем:

1. Генерируем начальное случайное решение (или получаем feasible при помощи эвристик);
2. Задаем начальную «температуру» – некий глобальный метапараметр;
3. Выполняем отжиг заданное число итераций;
4. Выполняем случайную модификацию решения;
5. Если значение функции стоимости для нового решения лучше, чем для старого, то принимаем его;
6. Если нет, то все равно можем принять новое решение, но лишь с некоторой вероятностью, которая тем больше, чем выше температура и чем ближе друг к другу по оценке старое и новое решение.

В качестве функции, которая дает нам вероятность принять/отклонить новое решение, будем использовать распределение Больцмана:

Видно, что эта величина может быть больше единицы в случае, когда новое решение лучше старого, но для нас это не проблема – это просто будет значить, что мы точно принимаем новое решение!

Эффективность метода повышается, если функция f(x) изменяется плавно при небольших изменениях x. Однако важно понимать, что выбор конкретных эвристик не универсален: все составляющие алгоритма — начальное состояние, определение «соседа», закон охлаждения — сильно зависят от характера задачи и подбираются экспериментально.

Пример: задача о N ферзях

Дана шахматная доска N × N и N ферзей. Нужно расставить их так, чтобы они не били друг друга.

Переведём задачу на язык отжига:

• Состояние x — список из N чисел, где x[i] — номер строки ферзя в столбце i. По одному ферзю в каждом столбце — вертикальные бои исключены самим представлением.
• Целевая функция f(x) — число пар ферзей, бьющих друг друга (по строкам и диагоналям). Решению соответствует f(x) = 0 — глобальный минимум известен заранее, что удобно для проверки.
• Соседнее состояние — один случайный ферзь переставляется в случайную новую строку своего столбца.

import random
import math

N = 8

def conflicts(queens):
    """Число пар ферзей, бьющих друг друга."""
    cnt = 0
    for i in range(N):
        for j in range(i + 1, N):
            if queens[i] == queens[j]:                # одна строка
                cnt += 1
            elif abs(queens[i] - queens[j]) == j - i:  # одна диагональ
                cnt += 1
    return cnt

def get_new_state(queens):
    """Сосед: один случайный ферзь — в случайную новую строку."""
    new_q = queens.copy()
    col = random.randrange(N)
    row = random.randrange(N)
    while row == new_q[col]:
        row = random.randrange(N)
    new_q[col] = row
    return new_q

def simulated_annealing(t_0=100.0, t_min=0.005):
    queens = [random.randrange(N) for _ in range(N)]
    cost = conflicts(queens)
    t_k = t_0
    k = 1
    while t_k > t_min and cost > 0:
        t_k = t_0 / (1 + k)      # плавное охлаждение температуры
        k += 1
        new_q = get_new_state(queens)
        new_cost = conflicts(new_q)
        dcost = new_cost - cost
        # лучше — принимаем всегда; хуже — с вероятностью exp(-dcost / t_k)
        if dcost <= 0 or random.random() < math.exp(-dcost / t_k):
            queens, cost = new_q, new_cost
    return queens, cost, k

queens, cost, k = simulated_annealing()
print("Конфликтов:", cost, "| итераций:", k)
print("Расстановка:", queens)

#include <iostream>
#include <vector>
#include <cmath>
#include <random>

const int N = 8;

std::mt19937 gen(std::random_device{}());
std::uniform_real_distribution<> uni(0.0, 1.0);

// Число пар ферзей, бьющих друг друга
int conflicts(const std::vector<int>& q) {
    int cnt = 0;
    for (int i = 0; i < N; ++i)
        for (int j = i + 1; j < N; ++j) {
            if (q[i] == q[j]) ++cnt;                        // одна строка
            else if (std::abs(q[i] - q[j]) == j - i) ++cnt; // одна диагональ
        }
    return cnt;
}

// Сосед: один случайный ферзь — в случайную новую строку
std::vector<int> get_new_state(const std::vector<int>& q) {
    std::vector<int> nq = q;
    int col = gen() % N;
    int row = gen() % N;
    while (row == nq[col]) row = gen() % N;
    nq[col] = row;
    return nq;
}

int main() {
    std::vector<int> queens(N);
    for (int i = 0; i < N; ++i) queens[i] = gen() % N;
    int cost = conflicts(queens);

    double t_0 = 100.0, t_min = 0.005, t_k = t_0;
    int k = 1;
    while (t_k > t_min && cost > 0) {
        t_k = t_0 / (1 + k);          // плавное охлаждение
        ++k;
        std::vector<int> nq = get_new_state(queens);
        int ncost = conflicts(nq);
        int dcost = ncost - cost;
        if (dcost <= 0 || uni(gen) < std::exp(-dcost / t_k)) {
            queens = nq;
            cost = ncost;
        }
    }

    std::cout << "Конфликтов: " << cost << " | итераций: " << k << "\n";
    std::cout << "Расстановка:";
    for (int r : queens) std::cout << " " << r;
    std::cout << "\n";
    return 0;
}

Алгоритм находит корректную расстановку в подавляющем большинстве запусков — как правило, за несколько сотен итераций (сравните с полным перебором: даже при одном ферзе в каждом столбце вариантов 8⁸ ≈ 16,8 миллиона). Изредка поиск «застревает» с одним-двумя конфликтами — температура успевает упасть раньше, чем найден выход из локального минимума; в таком случае алгоритм просто перезапускают. Это и есть плата за отсутствие гарантий, о котором говорилось в определении эвристического алгоритма.

3 Муравиный алгоритм

Муравьиный алгоритм (Ant Colony Optimization, ACO) — это метаэвристический метод для решения задач оптимизации, вдохновленный поведением муравьев в поиске пищи. Алгоритм был разработан в 1992 году Марко Дориги и его коллегами, и с тех пор широко применяется для решения различных задач, таких как задача коммивояжера (TSP), расписания, оптимизация маршрутов и др.

Основные идеи

1. Поведение муравьев: Муравьи находит кратчайший путь к пище, следуя за феромонами, которые они оставляют на пути.
2. Феромоны: Муравьи оставляют следы феромонов на пройденном пути. Эти следы испаряются со временем, но феромоны остаются на пути, по которому прошел муравей.
3. Алгоритм: Исключает полное переборное решение. Вместо этого множество агентов (муравьев) исследуют пространство решений параллельно, используя феромоны для принятия решений.
4. Использование феромонов: Пути, по которым проходит больше муравьев, становятся более привлекательными для других муравьев, так как оставшиеся феромоны усиливаются.

Шаги работы алгоритма

1. Инициализация: Инициализация феромонов на всех путях с равными значениями.
2. Построение решения: Каждый муравей строит решение, выбирая следующий шаг на основе вероятности, которая зависит от количества феромонов и длины пути.
3. Обновление феромонов:
   • Феромоны испаряются со временем (уменьшаются).
   • Добавляются новые феромоны на пути, пройденные муравьями. Чем лучше путь, тем больше феромонов добавляется.
4. Повторение: Алгоритм повторяется, пока не будет найдено оптимальное решение или не исчерпано количество итераций.

Параметры алгоритма

1. Количество муравьев: Количество агентов, которые будут искать решение.
2. Испарение феромонов: Коэффициент, определяющий скорость, с которой феромоны исчезают со времени.
3. Параметры α и β: Они контролируют влияние феромонов (α) и расстояния (β) на выбор пути.
4. Число итераций: Количество циклов, в которых муравьи будут искать решение.

Муравьиный алгоритм применим для различных задач:

• Задача коммивояжера (TSP): Поиск кратчайшего пути, который посещает все города.
• Оптимизация маршрута: Для логистики и доставки товаров.
• Распределение задач: Например, при решении задачи о распорядке работы оборудования.

Пример использования

Пример использования муравьиного алгоритма для решения задачи коммивояжера (TSP) приведен ниже.

import random
import numpy as np

# Матрица расстояний (граф с 4 вершинами)
distances = np.array([
    [0, 2, 2, 5],
    [2, 0, 3, 4],
    [2, 3, 0, 1],
    [5, 4, 1, 0]
])

num_ants = 5
num_iterations = 100
alpha = 1      # Влияние феромона
beta = 2       # Влияние расстояния
evaporation = 0.5
pheromone_deposit = 100

n = len(distances)
pheromones = np.ones((n, n))  # начальный феромон

def choose_next_city(current, visited):
    probabilities = []
    for i in range(n):
        if i in visited:
            probabilities.append(0)
        else:
            prob = (pheromones[current][i] ** alpha) * ((1.0 / distances[current][i]) ** beta)
            probabilities.append(prob)
    total = sum(probabilities)
    probabilities = [p / total for p in probabilities]
    return random.choices(range(n), weights=probabilities)[0]

def run_ant():
    path = [random.randint(0, n - 1)]
    while len(path) < n:
        next_city = choose_next_city(path[-1], path)
        path.append(next_city)
    return path

def path_length(path):
    return sum(distances[path[i]][path[i + 1]] for i in range(len(path) - 1)) + distances[path[-1]][path[0]]

for iteration in range(num_iterations):
    all_paths = [run_ant() for _ in range(num_ants)]
    all_lengths = [path_length(p) for p in all_paths]

    # Обновляем феромоны
    pheromones *= (1 - evaporation)
    for i in range(num_ants):
        path = all_paths[i]
        length = all_lengths[i]
        for j in range(len(path) - 1):
            a, b = path[j], path[j + 1]
            pheromones[a][b] += pheromone_deposit / length
            pheromones[b][a] += pheromone_deposit / length

# Выводим лучший путь
best_path = min([run_ant() for _ in range(100)], key=path_length)
print("Лучший путь:", best_path)
print("Длина:", path_length(best_path))

#include <iostream>
#include <vector>
#include <cmath>
#include <limits>
#include <random>
#include <algorithm>

using namespace std;

const int N = 4; // количество городов
const int NUM_ANTS = 5;
const int NUM_ITERATIONS = 100;
const double ALPHA = 1.0;     // Влияние феромона
const double BETA = 2.0;      // Влияние расстояния
const double EVAPORATION = 0.5;
const double Q = 100.0;

double distances[N][N] = {
    {0, 2, 2, 5},
    {2, 0, 3, 4},
    {2, 3, 0, 1},
    {5, 4, 1, 0}
};

double pheromones[N][N];

random_device rd;
mt19937 gen(rd());

int select_next_city(int current, const vector<bool>& visited) {
    vector<double> probabilities(N, 0.0);
    double sum = 0.0;

    for (int i = 0; i < N; ++i) {
        if (!visited[i]) {
            double tau = pow(pheromones[current][i], ALPHA);
            double eta = pow(1.0 / distances[current][i], BETA);
            probabilities[i] = tau * eta;
            sum += probabilities[i];
        }
    }

    uniform_real_distribution<> dis(0.0, sum);
    double r = dis(gen);
    double cumulative = 0.0;

    for (int i = 0; i < N; ++i) {
        if (!visited[i]) {
            cumulative += probabilities[i];
            if (r <= cumulative) return i;
        }
    }

    // fallback (должно не понадобиться)
    for (int i = 0; i < N; ++i)
        if (!visited[i]) return i;

    return 0;
}

vector<int> run_ant() {
    uniform_int_distribution<> dis(0, N - 1);
    int start = dis(gen);
    vector<int> path = {start};
    vector<bool> visited(N, false);
    visited[start] = true;

    while (path.size() < N) {
        int next = select_next_city(path.back(), visited);
        path.push_back(next);
        visited[next] = true;
    }

    return path;
}

double path_length(const vector<int>& path) {
    double length = 0.0;
    for (int i = 0; i < path.size() - 1; ++i)
        length += distances[path[i]][path[i + 1]];
    length += distances[path.back()][path[0]]; // возвращение в начало
    return length;
}

int main() {
    // Инициализация феромонов
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            pheromones[i][j] = 1.0;

    for (int iter = 0; iter < NUM_ITERATIONS; ++iter) {
        vector<vector<int>> all_paths(NUM_ANTS);
        vector<double> lengths(NUM_ANTS);

        for (int k = 0; k < NUM_ANTS; ++k) {
            all_paths[k] = run_ant();
            lengths[k] = path_length(all_paths[k]);
        }

        // Испарение феромонов
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j)
                pheromones[i][j] *= (1.0 - EVAPORATION);

        // Добавление феромонов
        for (int k = 0; k < NUM_ANTS; ++k) {
            double contribution = Q / lengths[k];
            const auto& path = all_paths[k];
            for (int i = 0; i < N - 1; ++i) {
                int a = path[i];
                int b = path[i + 1];
                pheromones[a][b] += contribution;
                pheromones[b][a] += contribution;
            }
            // возврат в начало
            pheromones[path.back()][path[0]] += contribution;
            pheromones[path[0]][path.back()] += contribution;
        }
    }

    // Поиск лучшего пути после обучения
    double best_length = numeric_limits<double>::max();
    vector<int> best_path;

    for (int i = 0; i < 100; ++i) {
        auto path = run_ant();
        double len = path_length(path);
        if (len < best_length) {
            best_length = len;
            best_path = path;
        }
    }

    cout << "Лучший путь: ";
    for (int city : best_path)
        cout << city << " ";
    cout << best_path[0] << endl; // возвращение
    cout << "Длина: " << best_length << endl;

    return 0;
}

Для матрицы из листинга оптимальный замкнутый маршрут имеет длину 9 (это легко проверить полным перебором — всего 3! = 6 существенно различных маршрутов), и муравьиный алгоритм устойчиво его находит.

4 Вопросы для закрепления

5 Список использованных источников и благодарностей

1. Алгоритмика
2. QMLCourse