Мне нравится сайт Хабр, особенно нравится как люди умеют грамотно излагать свои «технические» мысли на бумаге, чего не могу сказать о себе. Поэтому максимально постарался все описать в комментариях — мне так намного легче, да и Вам удобнее будет разобраться.
Код написан на MATLAB. При определенном написании синтаксис понять несложно. При написании алгоритмов старался сохранить баланс между «пониманием» кода и скоростью выполнения, исключив векторизацию, а где она есть — в комментариях написан альтернативный код. Если возникнут вопросы — постараюсь ответить в комментариях. Это первая часть серий «Оптимизация на примере», дальнейшие публикации будут выходить в зависимости от Ваших комментариев. То есть если нужно сравнение с другими эвристическими методами, либо если у Вас возникнет какое-то предложение (скажем изменить формулу температуры, либо начальное расположение муравьев), где можно доработать, а где изменить часть алгоритма — пишите. Будет интересно всем, да и мне тоже. Я же по возможности реализую и также опубликую.
Цель публикации, в первую очередь, показать на подробных комментариях в коде принцип работы двух алгоритмов, сравнение — второстепенная цель. В данной статье сравнить «чисто» не получится, так как в двух алгоритмах, особенно в муравьином, содержится приличное число начальных изменяемых параметров. Более того, возникли некоторые вопросы, но об этом чуть позже. В имитационном отжиге проще, однако, тоже есть свои нюансы. Для тех, кто только знакомится с данными алгоритмами рекомендую сначала ознакомиться с полезным материалом, расположенным в конце статьи.
Теория двух алгоритмов есть на Хабре да и на других сайтах, однако открытого кода муравьиного алгоритма с поясняющими комментариями толком не нашел, есть псевдокод, есть отдельные его части. Код же имитационного отжига на MATLAB был выложен ранее, однако, для сравнения с другими видами оптимизации у него небольшая скорость. Поэтому написал свой, увеличив скорость в разы. Оба алгоритма написал без встроенных собственных функций, потому что при небольшом коде, по-моему мнению, так воспринимать структуру алгоритма легче. В местах, где можно запутаться старался дать более развернутые комментарии. Исходя из этого, в теорию вдаваться не буду, приведу лишь следующее:
Имитационный отжиг
Метод имитации отжига отражает поведение материального тела при управляемом охлаждении. Как показали исследования, при отвердевании расплавленного материала его температура должна уменьшаться постепенно, вплоть до момента полной кристаллизации. Если процесс остывания протекает слишком быстро, образуются значительные нерегулярности структуры материала, которые вызывают внутренние напряжения. Быстрая фиксация энергетического состояния тела на уровне выше нормального аналогична сходимости оптимизационного алгоритма к точке локального минимума. Энергия состояния тела соответствует целевой функции, а абсолютный минимум этой энергии – глобальному минимуму [1].
Муравьиный алгоритм
Один из эффективных полиномиальных алгоритмов для нахождения приближённых решений задачи коммивояжёра, а также решения аналогичных задач поиска маршрутов на графах. Суть подхода заключается в анализе и использовании модели поведения муравьёв, ищущих пути от колонии к источнику питания и представляет собой метаэвристическую оптимизацию [2].
Есть множество модификаций муравьиных алгоритмов. Где-то разное обновление феромонов, где-то имеются «элитные» муравьи, а где-то различные расположения муравьев при старте. В коде просто что-то закомментируйте, что-то раскомментируйте.
% метод отжига ( на примере задачи коммивояжера )
%--------------------------------------------------------------------------
tic
% clearvars -except cities
clearvars
% -----------------------------ИТЕРАЦИИ------------------------------------
% кол-во итераций
m = 100000;
% -----------------------------ПАРАМЕТРЫ-----------------------------------
% начальная температура
Tstart = 10000;
% конечная температура
Tend = 0.1;
% начальная температура для вычислений
T = Tstart;
% расстояние
S = inf;
% количество городов
n = 100;
% --------------------------------ПАМЯТЬ-----------------------------------
% матрица расстояний
dist = zeros(n,n);
% -------------------------------------------------------------------------
% генерация городов (x,y)
cities = rand(n,2)*100;
% задаем случайный маршрут
ROUTE = randperm(n);
% создаем матрицу расстояний
for i = 1:n
for j = 1:n
% dist ( расстояния )
dist(i,j) = sqrt((cities(i,1) - cities(j,1))^2 + ...
(cities(i,2) - cities(j,2))^2);
end
end
% поехали оптимизировать, время от кол-ва итераций
for k = 1:m
% сбрасываем потенциальное расстояние
Sp = 0;
% здесь условие создания потенциальных маршрутов, ROUTEp -
% потенциальный маршрут
% потенциальный маршрут
ROUTEp = ROUTE;
% два случайных города
transp = randperm(n,2);
% переворот вектора, взяли два сгенерированных числа transp и
% перевернули маршрут между ними
if transp(1) < transp(2)
ROUTEp(transp(1):transp(2)) = fliplr(ROUTEp(transp(1):transp(2)));
else
ROUTEp(transp(2):transp(1)) = fliplr(ROUTEp(transp(2):transp(1)));
end
% вычисляем энергию (расстояние) потенциального маршрута
for i = 1:n-1
Sp = Sp + dist(ROUTEp(i),ROUTEp(i+1));
end
Sp = Sp + dist(ROUTEp(1),ROUTEp(end));
% если он меньше то, потенциальный маршрут становится основным ...
% если нет, смотрим, осуществляется ли переход
if Sp < S
S = Sp;
ROUTE = ROUTEp;
else
% случайно выбранное число от 0 до 1
RANDONE = rand(1);
% вычисляем вероятность перехода
P = exp((-(Sp - S)) / T);
if RANDONE <= P
S = Sp;
ROUTE = ROUTEp;
end
end
% уменьшаем температуру
T = Tstart / k;
% проверяем условие выхода
if T < Tend
break;
end;
end
% рисуем графику
citiesOP(:,[1,2]) = cities(ROUTE(:),[1,2]);
plot([citiesOP(:,1);citiesOP(1,1)],[citiesOP(:,2);citiesOP(1,2)],'-r.')
% очищаем переменые
clearvars -except cities ROUTE S
% смотрим время
toc
% муравьиный алгоритм ( на примере задачи коммивояжера )
% -------------------------------------------------------------------------
tic
% clearvars -except cities
clearvars
% ------------------------------ИТЕРАЦИИ-----------------------------------
% кол-во итераций ( поколений )
age = 50;
% кол-во муравьев в поколении
countage = 30;
% кол-во городов
n = 30;
% ------------------------------ПАРМЕТРЫ-----------------------------------
% альфа - коэффициент запаха, при 0 будем ориентироваться только на
% кратчайший путь
a = 1;
% бета - коэффициент расстояния, при 0 будем
% ориентироваться только на оставляемый запах
b = 2.5;
% коэффициент испарения ( память поколений )
e = 0.5;
% количество выпускаемых феромонов
Q = 1;
% кол-во элитных муравьев ( увеличение лучшей дистанции в elite раз )
elite = 1;
% начальный феромон
ph = 0.01;
% -------------------------------ПАМЯТЬ------------------------------------
% матрица расстояний
dist = zeros(n,n);
% матрица обратных расстояний
returndist = zeros(n,n);
% матрица маршрута муравьев в одном поколении
ROUTEant = zeros(countage,n);
% вектор расстояний муравьев в одном поколении
DISTant = zeros(countage,1);
% вектор лучших дистанций на каждой итерации
bestDistVec = zeros(age,1);
% лучший начальный маршрут
bestDIST = inf;
% оптимальные маршруты
ROUTE = zeros(1,n+1);
% перестановка городов без повторений ( для выхода муравьев )
RANDperm = randperm(n);
% матрица вероятностей
P = zeros(1,n);
% -------------------------------------------------------------------------
% генерация городов (x,y)
cities = rand(n,2)*100;
% матрица начальных феромонов
tao = ph*(ones(n,n));
% создаем матрицу расстояний и матрицу обратных расстояний
for i = 1:n
for j = 1:n
% dist ( расстояния )
dist(i,j) = sqrt((cities(i,1) - cities(j,1))^2 + ...
(cities(i,2) - cities(j,2))^2);
% nn ( обратные расстояния )
if i ~= j
returndist(i,j) = 1/sqrt((cities(i,1) - cities(j,1))^2 + ...
(cities(i,2) - cities(j,2))^2);
end
end
end
% итерации
for iterration = 1:age
% муравьи ( одно поколение)
for k = 1:countage
% ****************** НАЧАЛЬНОЕ РАСПОЛОЖЕНИЕ МУРАВЬЕВ ******************
% выбирайте какой нужно
% каждый муравей располагается случайно
% ROUTEant(k,1) = randi([1 n]);
% с каждого города выходит один муравей ( без совпадений ), кол-во
% городов и кол-во муравьев в поколении должны быть равны
ROUTEant(k,1) = RANDperm(k);
% с конкретного города выходят все муравьи в данном случа с 1-ого
% ROUTEant(k,1) = 1;
% тут маршрут первому поколению задаем либо произвольный, либо с каждого
% города разный, либо с одного города все, а следующее поколение выходит по
% концам первых
% if iterration == 1
% ROUTEant(k,1) = randi([1 n]);
% % ROUTEant(k,1) = RANDperm(k);
% % ROUTEant(k,1) = 1;
% else
% ROUTEant(k,1) = lastROUTEant(k);
% end
% *********************************************************************
% путь каждого муравья, начиная со второго, так как первый выбран
for s = 2:n
% полуаем индекс выбранного города
ir = ROUTEant(k,s-1);
% вероятность посещения городов ( числитель ) , в числителе у нас
% следующее: tao^a*(1/S)^b
% 1/S -это returndist.
% поскольку данное значение будет повторяться (кол-во муравьев * на
% колонию * кол-во городов) раз, то еще один цикл писать не выгодно,
% скорость работы при таких вычислениях падает. Поэтому написал в
% этом моменте векторно. На обычном языке будет так:
% for c = 1:n
% P(1,c) = tao(ir,c).^a * returndist(ir,c).^b;
% end
P = tao(ir,:).^a .* returndist(ir,:).^b;
% получили числители (в формуле вероятности перехода к k-ому городу)
% для n городов, однако в некоторых мы уже побывали, нужно исключить
% их
% проставляем нули в числитель туда, где уже были, чтобы
% вероятность перехода была 0, следовательно в сумме знаменателя
% формулы данный город учитываться не будет
P(ROUTEant(k,1:s-1)) = 0;
% получаем вероятности перехода ( сумма строк должна быть = 1 )
P = P ./ sum(P);
% смотрим в какой город осуществляется переход
RANDONE = rand;
getcity = find(cumsum(P) >= RANDONE, 1, 'first');
% если с вероятностями что-то не так, выдача ошибки
if isempty(getcity) == 1
disp 'Ошибка в вероятностях!'
end
% присваем s-ый город в путь k-ому муравью
ROUTEant(k,s) = getcity;
end
% получаем маршрут k-ого муравья
ROUTE = [ROUTEant(k,1:end),ROUTEant(k,1)];
% сброс длины
S = 0;
% вычисляем маршрут k-ого муравья
for i = 1:n
S = S + dist(ROUTE(i),ROUTE(i+1));
end
% путь k-ого муравья, массив дистанций k-ых муравьев
DISTant(k) = S;
% присваевыем лучший маршрут и S
if DISTant(k) < bestDIST
bestDIST = DISTant(k);
bestROUTE = ROUTEant(k,1:end);
end
% вектор "последних" городов k-ых муравьев ( выбирается для старта
% муравьев нового поколения с тех городов, где закончили путь
% предыдущее поколение)
% lastROUTEant = ROUTEant(1:end,end);
end
% -----------------------ИЗМЕНЯЕМ КОЛИЧЕСТВО ФЕРОМОНОВ---------------------
% Испаряем феромоны "старого" пути е - коэффициент испарения
tao = (1-e)*tao;
for u = 1:countage
% получаем маршрут u-ого муравья
ROUTE = [ROUTEant(u,1:end), ROUTEant(u,1)];
% для каждого города
for t = 1:n
x = ROUTE(t);
y = ROUTE(t+1);
% считаем новый феромон
if DISTant(u) ~= max(DISTant)
% тут по формуле добавки феромона в зависимости от
% дистанции
tao(x,y) = tao(x,y) + Q / DISTant(u);
else
% усиливаем ребера лучшего пути
tao(x,y) = tao(x,y) + (elite*Q) / DISTant(u);
end
end
end
end
% строим графику
citiesOP(:,[1,2]) = cities(bestROUTE(:),[1,2]);
plot([citiesOP(:,1);citiesOP(1,1)],[citiesOP(:,2);citiesOP(1,2)],'.r-')
clearvars -except cities ROUTE bestDIST bestROUTE ROUTEant
toc
Итак, давайте сравним, сгенерируем для начала 30 городов.
Имитационный отжиг:
Параметры:
- Количество городов — 30
- Начальная температура — 10 000
- Конечная температура — 0,01
- Формула температуры — начальную температуру / k-ую итерацию
- Число итераций — 100 000
- Функция вероятности принятия — exp(-?E/T)
- Определение потенциального маршрута (порождающего семейства) — разворот части вектора (текущего маршрута) от двух случайно выбранных чисел равномерным распределением
Теперь посмотрим на работу муравьиного алгоритма на тех же координатах.
Муравьиный алгоритм:
Параметры:
- Количество итераций (поколений) — 60
- Количество муравьев в поколении — 30
- Количество городов — 30
- альфа (коэф. ориентации на феромоны) — 0,85
- бета (коэф. ориентации на длину пути) — 3,8
- е (коэф. испарение феромонов) — 0,5
- elite (количество элитных муравьев) — нет
- расположение муравьев в первых городах — разные, без повторений. Кол-во городов = числу муравьев в поколении
- обновление феромонов — глобальное (феромоны обновляются после прохода 1-ого поколения)
- вероятность перехода из одной вершины в другую —
Перед Вами наглядная работа двух алгоритмов. Еще раз отмечу, что подробное объяснение (особенно муравьиного алгоритма) в комментариях кода.
Обратите внимание, что в имитационном отжиге оптимум нашелся уже на 20000 итерации, дальше, когда температура стала достаточно низкой, вероятность перехода на новые маршруты значительно снизилась, поэтому алгоритм «встал». Тут можно ввести критерии остановки. Чем выше начальная температура, тем больше вероятность попасть на оптимальный маршрут, и тем больше итераций будет совершено, однако, далеко не факт, что мы найдем глобальный оптимум.
В муравьином алгоритме же все иначе. Главный плюс муравьиного алгоритма заключается в том, что он никогда не остановится, он всегда будет искать лучшие пути. При бесконечном числе итераций (поколений) вероятность найти глобальный оптимум будет стремиться к 1. Однако, уже в самом начале встал вопрос: какие начальные параметры выбрать? Ответ: методом подбора. Снова оптимизируем. Данные параметры муравьиной оптимизации оптимизировал с помощью генетического алгоритма не на данных координатах. Там были другие координаты и количество городов было 50. Толком не получив хорошей сходимости, получил некоторые параметры, которые примерно соответствуют на данном примере.
На одном тесте далеко не уехать, поэтому привожу таблицу из серий тестов.
Видно стабильность двух алгоритмов, однако разброс от среднего значения у муравьиного алгоритма меньше.
Теперь посмотрим работу двух алгоритмов на 100 сгенерированных городов.
Имитационный отжиг:
Параметры:
- Количество городов — 100
- Начальная температура — 42 000
- Конечная температура — 0,01
- Формула температуры — начальную температуру / k-ую итерацию
- Число итераций — 420 000
- Функция вероятности принятия — exp(-?E/T)
- Определение потенциального маршрута (порождающего семейства) — разворот части вектора (текущего маршрута) от двух случайно выбранных чисел равномерным распределением
Теперь посмотрим на работу муравьиного алгоритма на тех же координатах.
Муравьиный алгоритм:
Параметры:
- Количество итераций (поколений) — 50
- Количество муравьев в поколении — 100
- Количество городов — 100
- альфа (коэф. ориентации на феромоны) — 1
- бета (коэф. ориентации на длину пути) — 3
- е (коэф. испарение феромонов) — 0,5
- elite (количество элитных муравьев) — нет
- расположение муравьев в первых городах — разные, без повторений. Кол-во городов = числу муравьев в поколении
- обновление феромонов — глобальное (феромоны обновляются после прохода 1-ого поколения)
Здесь параметры взял «на глаз».
Снова муравьиный алгоритм показал результат получше, однако, давайте снова посмотрим на таблицу из серий тестов:
Здесь далеко ушел вперед имитационный отжиг. Вообще, до того как я прочитал статью, к имитационному отжигу относился скептически. В разработке торговых систем (в чем я замешан) и их оптимизации метод отжига не подходит, так как нам важен не сам локальный оптимум, а «разброс» вокруг него. На условном примере: есть скользящая средняя в 15 дней на какой-то бумаге, за год доходность показала 10000, а при 14 днях (-1000). Здесь нам важно, чтобы около параметра оптимального равномерно убывала прибыль. Поэтому в разных ситуациях — разная оптимизация.
Что касается задачи коммивояжера, то я, пусть даже и не в «чистом» соревновании (потому как правильно подобрать параметры для муравьиного алгоритма не вышло), отдаю победу имитационному отжигу. При добавлении критерия остановки, тех же расстояний без особого изменения стандартного отклонения можно добиться за 2,5 секунды, чего не скажешь о муравьином алгоритме. По серии тестов, значительная часть которых не попала сюда, муравьиный алгоритм блестяще работает до 50 городов, дальше начинает «моросить», время оптимизации сильно замедляется, что в принципе и логично. Однако он уверено идет к глобальному оптимуму, пусть даже и медленно.
Где-то в глубине души все равно хочется, чтобы муравьиный алгоритм проявил себя лучше. Может быть поэтому не покидает чувство, что допустил где-то ошибку, и не одну. Поэтому те, кого заинтересовал, и те, кто только начинает путь в этом направлении — присоединяйтесь. Пишите комментарии с предложениями, как по поводу алгоритмов, так и по поводу того, как лучше подавать материал.
А пока посмотрим на сегодняшнего победителя. 500 городов, 86 секунд оптимизации (и это далеко не предел).
Очень круто!
Догонит или нет муравьиный алгоритм? Обгонит ли его генетический? Кто будет эффективнее? Об этом и другом в следующих публикациях.
Полезный материал:
» Видео про муравьиный алгоритм (на русском)
» Видео про имитационный отжиг (на русском)
» Статья про муравьиный алгоритм (на русском)
» Статья про имитационный отжиг ( на русском)
Комментарии (14)
QtRoS
31.08.2016 23:17О, с удовольствием вспоминаю курс Discrete Optimization от Pascal van Hentenryck. Там в последнем задании были тысячи городов, и если кому-то интересно, могу выложить картинку, которая получилась. Одно из популярных решений (кстати, российского участника курсов) даже запечатлили на майке:
daiver19
01.09.2016 00:07+3Мини-поправка: Имитация отжига, а не «имитационный отжиг».
Отжиг — это, пожалуй, самый универсальный и быстрый алгоритм оптимизации. Очень большой процент TopCoder-овских марафонов был выигран отжигом, да я и сам использовал его в куче задач. Это может удивлять поначалу, но стоит это принять и стараться использовать в подходящих ситуациях.Grossmend
01.09.2016 00:46Если бы метода отжига не существовало и мне бы предложили написать данный алгоритм, то я бы даже не брался. Как-то сложно поверить в него. Однако, он выдает отличные результаты. Вы правильно подметили что он универсален.
daiver19
01.09.2016 00:49+1Да, я понимаю. Чтобы прочувствовать это стоит попробовать простой hill climbing (случайно меняем решение и если оно улучшилось, то оставляем его), который показывает весьма неплохие результаты. Отжиг кажется логичным развитием этого алгоритма.
Undiabler
01.09.2016 00:21Лично мне интуитивно кажется что ни первый ни второй алгоритм не работают оптимально.
Задача коммивояжера это оптимальный путь, значит, при наличие «скоплений» городов в определенных зонах, алгоритм должен в первую очередь искать оптимальные пути выбирая только из близлежащих городов, и не учитывать удаленные точки.
Что в одном что во втором примере проверяются маловероятные связки и это кажется мне не оптимальным. Поправьте если ошибаюсь.Grossmend
01.09.2016 00:41В Ваших словах есть истина. Но это будет совсем другой алгоритм, здесь же сравниваются только эти два. Есть идея сделать как Вы написали. Сначала найти «скопления», там найти локальный оптимум, затем уже найти оптимум между «скоплениями». Если заметить, то как раз таки в муравьином алгоритме, уже с первых итераций в «пучках» образуются маршруты, затем алгоритм «танцует» около них.
zartdinov
01.09.2016 00:21+1Не специалист, но вроде можно сильно ускорить вашу реализацию метода отжига, если в каждой итерации не пересчитывать расстояние маршрута полностью. Если правильно представляю, то достаточно сравнить длину двух старых и двух новых отрезков чтобы понять уменьшится общая протяженность или нет. В случае большого количества городов можно потратить это время на новые попытки улучшить результат.
Муравьиный алгоритм еще не смотрел. Спасибо за статью, очень интересно.Grossmend
01.09.2016 00:42Здравствуйте, можно чуть подробней.
daiver19
01.09.2016 00:57+1Когда мы меняем две вершины местами остальные вершины остаются на прежнем месте. Поэтому операция перемены мест двух вершин эквивалентна удалению 3-4 ребер из пути и добавлению новых 3-4 ребер. Значит мы можем пересчитывать только эти ребра, итого модификация занимает константное время. Собственно этот подход используют как раз из-за скорости вычисления. Если пересчитывать весь маршрут, то можно делать очень сложные модификации, но обычно это не нужно. Хотя добавление более сложных операций вроде вставки вершины из одного места в другое может несколько ускорить сходимость (ценой усложнения реализации эффективного рассчета пути).
Asker1987
01.09.2016 14:47Отдавать «победу имитационному отжигу» (ИО) довольно смелый поступок, особенно, в задаче коммивояжера. Либо нужно оговорить, в каких условиях ИО выигрышна. При серьезных временных ограничениях многие алгоритмы «побеждают» муравьиный алгоритм(МА), но чаще важно найти именно оптимум. Можно говорить о конкретных реализациях алгоритмов, но если хотите сравнить, то возьмите научную литературу, статьи уважаемых зарубежных авторов, в открытых источниках приведены бенчмарки. На сколько мне известно, то в этой задаче правят балом генетический алгоритма и МА.
По статье и реализации МА:
1. Нет ссылок на работы других авторов. Конечно, эта статья не для научного журнала, но если делаете громкие заявления, то нужно делать и громкие исследования. Ведь МА и ИО не вчера разработаны, их давно сравнили, в т.ч. на задаче коммивояжера.
2. У Вас приведен пример простейшего МА. Нет смысла говорить о «победе» одного научного направления над другим. Существуют десятки, а то и больше модификаций МА. Простой МА лишь отражает основную идею использования «глобальной памяти» через феромоны.
3. С числом муравьев допущена типовая ошибка — в общем случае число муравьев не равно числу точек. Размер колонии может быть любой и муравьи рандомно располагаются. Иначе МА «моросит» при большом кол-ве точек. Число муравьев должно коррелировать с числом ребер(!), а не вершин. Сравним с естественной средой: чем больше муравьев, тем больше площадь они могут исслдеовать, НО площадь исследования в графе — это число ребер, а не вершин. Слишком вероятна стагнация.
4. Коэффициент бета очень не желательно делать дробной. Честно, не знаю, как устроен матлаб, в С++ всегда избегал дробного возведения в степень, т.к. это увеличивает время в МА, порой даже в разы.
PS. Код не смотрел, уверен там все филигранно)
Grossmend
01.09.2016 22:11Победу отдал на «сегодняшний день». Пункты 2-3, а также другие модификации предназначаются для следующих публикаций. 1 – учту. Спасибо.
saboteur_kiev
> А пока посмотрим на сегодняшнего победителя. 500 городов, 86 секунд оптимизации (и это далеко не предел).
> Очень круто!
да!