Оптимизация на примере. Муравьиный алгоритм (ACS) против Метода отжига. Часть 2
Данная статья не доказывает какой метод лучше в общем, так как у обоих алгоритмов есть множество модификаций, а также множество потенциальных модификаций. Здесь всего лишь определяется победитель на определенной модификации и на определенных параметрах. Цель же данной статьи — расширение кругозора читателей в области дискретной оптимизации, и, надеюсь, предложения по улучшению работы обоих алгоритмов.
Модификация муравьиного алгоритма выбрана ACS (Ant Colony System), которую предложили Марко Дориго и Лука Гамбарделла в 1997 году. Главные отличия от AS (Ant System) это:
1) задается баланс между самым привлекательным городом, и выбором как в обычном AS
arg max{[τ (r, u)]^α*[ω (r, u)]^β} если q <= q0 (Формула 1)
u ϵ Jk«r»
В обратном случае выбираем переход по AS
, где [τ (r, u)] — уровень феромона на ребре (r, u), [ω (r, u)] — вес обратный расстоянию на ребре (r, u), β — регулируемый параметр, чем он выше, тем алгоритм будет склонен выбирать город, имеющий меньшее расстояние, α — равна 1, q — случайно выбранное число, q0 — вероятность выбора того, что переход из одной вершины в другую будет идти по формуле 1, u — города, еще не посещенные
2) помимо глобального обновления феромонов происходит еще и локальное. Уровень феромонов меняется при прохождении каждого муравья на итерации (тут ближе к естественной среде обитания муравьев)
τ (r, s)=(1-p)*τ (r, s)+p*τ0 (Формула 2)
, где p — уровень локального обновления, τ0 = значение начального феромона, который вычисляется следующим образом: τ0 = (n*Lnn)^-1, где Lnn — приблизительное оптимальное значение, которое может быть получено другим методом оптимизации.
3) при глобальном обновлении феромона, добавление происходит только к ребрам, либо лучших со времени начала работы алгоритма (глобальных лучших), либо к ребрам лучших на итерации (локальных лучших). Я применил на ребра глобальных лучших.
τ (r, s)=(1-e)*τ (r, s)+e*(Lbest^-1) (Формула 3)
, где e — уровень глобального обновления, Lbest — лучшая длина маршрута (самый короткий), либо на k-ой итерации, либо глобальный лучший.
Данные модификации дали значительную производительность алгоритму. Метод отжига оставил таким же, как и в прошлой статье, за исключением увеличения скорости. Спасибо читателю zartdinov за простое и очень эффективное предложение по увеличению скорости.
Теперь будем сравнивать задачи на известных координатах, таких как Oliver30, Eli51, Eli101, на которых найдено лучшее решение. Также выведем приближенные формулы временной зависимости двух алгоритмов от количества городов, и, наконец, попытаемся определить победителя на сегодня с учетом всех этих факторов.
Начнем с задачи Oliver30 лучшее решение — 423,7406
Параметры ACS:
- Количество итераций (поколений) — 2500*
- Количество муравьев в поколении — 7*
- Количество городов — 30*
- альфа (коэф. ориентации на феромоны) — 1
- бета (коэф. ориентации на длину пути) — 2
- p (коэф. обновления феромонов локальное) — 0,09
- e (коэф. обновления феромонов глобальное) — 0,09
- q (коэф. выбора самого привлекательного города) — 0,9
- начальное расположение муравьев — случайный
* — изменяемые параметры от размерности задачи
Пару слов о количестве муравьев. Ошибкой будет полагать, что чем больше муравьев, тем лучше. При большом количестве муравьев в колонии производительность значительно падает. Во-первых, сильно увеличивается время работы алгоритма. Во-вторых, происходит переизбыток феромонов, что приводит к аналогии в естественной среде, именуемой кругом смерти муравьев. Таким образом, алгоритм застревает на локальном оптимуме.
Точного определения количества муравьев в колонии пока нет, есть приблизительное вычисление как в [1]. Тут же я старался найти точку оптимума, при которой увеличение количества муравьев не улучшало решение, а их уменьшение снижало результат. Также, если упрощенно, экспериментировал с одной колонией муравьев на определенных итерациях. Как только средний уровень феромона начинал повышаться — брал это за основу оптимального количества муравьев.
Итак результаты задачи Oliver30:
Пара графиков:
Четвертый график показывает, что алгоритм продолжает искать альтернативные пути, не останавливаясь на достигнутом. При увеличении числа муравьев до 50–100 разброс средних дистанций поколения сокращается в пределах 20–30, что приводит к застреванию.
Полный перебор вариантов для симметричной задачи коммивояжера составляет (n-1)!/2 или
4 420 880 996 869 850 977 271 808 000 000 для данной задачи
100% результат, отличная работа ACS
Посмотрим на метод отжига.
Параметры:
- Количество городов — 30*
- Начальная температура — 35 000*
- Конечная температура — 0,1
- Формула температуры — начальную температуру / k-ую итерацию
- Число итераций — 350 000*
- Функция вероятности принятия — exp (-ΔE/T)
- Определение потенциального маршрута (порождающего семейства) — разворот части вектора (текущего маршрута) от двух случайно выбранных чисел равномерным распределением
* — изменяемые параметры от размерности задачи
Результаты:
Как по времени, так и по качеству на 30 городах выигрывает с большим отрывом ACS. Тестировал два алгоритма не только на данной задаче, но и на других 30 — ACS безусловно побеждает простой метод отжига.
Теперь задача на Eil51 на 51 город, лучшее решение 426, 7000 итераций, кол-во муравьев 9
Пара графиков:
На 4-ом графике добавлена линия регрессии. Вообще, задача Eil51 во многих зарубежных исследованиях тестировалась на куда большем числе итераций. Может быть поэтому глобальный оптимум не был найден, честно говоря, также тестировал на больших итерациях и максимум что был найден это 427.
Пользуясь случаем давайте посмотрим «в реальном режиме» на изменение феромонов от числа итераций, мне данная картинка довольно понравилась.
Без глобального обновления и с добавлением общего коэффициента испаряемости как в AS, для большего эффекта.
Смотрим на метод отжига при 2 500 000 итераций
Довольно неплохой результат, но все же ACS еще впереди.
Теперь задача Eil101 на 101 город, лучшее расстояние 629, 9000 итераций, кол-во муравьев 11
Пара графиков:
Здесь 9000 итераций, довольно мало для 100-ой задачи, однако сравним данный результат с отжигом на том же временном интервале при 7 000 000 итераций.
Довольно стабильный результат, лучше чем показал ACS. Но ACS позволяет сделать поиск более управляемым, чем обычный отжиг (хотя в последнем как минимум можно было бы ввести критерий остановки, но об этом в следующих статьях). Поэтому сегодня до 100 вершин по всем параметрам выигрывает ACS. Более того, сильно ускорить метод отжига, скорее всего, больше не удастся, в то время как оптимизировать код ACS еще как возможно. (В данном случае код плохо оптимизирован).
Повторюсь, что главный плюс ACS (как и других модификаций муравьиного алгоритма), в поиске глобального оптимума, при бесконечном числе итераций вероятность нахождения глобального лучшего стремится к 1. Вопрос конечно же во времени. Быстро, но примерно, либо долго, но точно. Поэтому предлагаю построить зависимость времени (при одинаковых числах итерациях на каждом алгоритме) от количества городов.
Количество итераций ACS — 10 000, муравьев — 10;
Количество итераций SA — 7 000 000.
Ого! Видим, что метод отжига занимает константное время, вне зависимости от количества вершин. Построив регрессию, определяем временную сложность для ACS.
t = 0.0044939x^2 + 0.72097x + 3.8225 (Формула 4)
, где x — количество городов, t — время выполнения ACS
Если же два алгоритма до 100 вершин примерно шли вровень как по времени, так и по качеству (с небольшим опережением ACS), то совсем грубо можно предположить что на 1000 городах, в 10 000 итераций на ACS и в 7 000 000 итераций на SA результат должен быть схожем.
Проверим.
ACS 200 городов (случайные города), время — 311,54 с.
SA 200 городов (те же что и выше), время 103,19 с.
Запустил последовательно. Какова вероятность того, что оба покажут один результат? Интересный момент вышел, может и сотые доли совпали? Но этого уже не узнать ни Вам, ни мне)
Вообщем по серии тестов, и по 300 вершин выходит примерно одинаковый результат, с ростом же выше в плюс уходит метод отжига.
При временном ограничении и количеством вершин больше ста лучше проявляет себя простой метод отжига нежели ACS. Повторюсь, что именно ACS, а не MMAS, ACS local search, либо иная модификация.
Но иногда нужно найти именно глобальный оптимум, тут мало кто может потягаться с Муравьиными алгоритмами.
Теперь пару слов об ускорении метода имитации отжига.
Как подсказал читатель daiver19, то, конечно же, не стоит пересчитывать на каждой итерации маршрут.
Есть условный маршрут:
1 2 3 4 5 6 7 8 9
Случайно отобрали два числа, пусть будет (2,5)
Теперь достаточно посчитать расстояние (1,2) и (5,6) и посчитать расстояние (1,5) и (2,6)
Однако на асимметричной задаче так не пройдет.
Благодаря этому количество городов не влияет на время выполнения алгоритма, плюс убрал функции fliplr, которые занимали приличное количество времени. Также были сформированы заранее массивы случайных чисел. Все это в разы увеличило скорость алгоритма с предыдущей статьей.
Также одному из читателей интересен был бы результат при перестановке двух вершин, а не инвертировании пути между ними. Давайте посмотрим результат на 101 городах Eil101 на 1 000 000 итераций.
Инверсия пути значительно лучше.
Много чего еще хотелось бы показать и рассказать, но и так статья выходит достаточно длинной. В следующих публикациях «наворотим» метод отжига, попытаемся сделать его более управляемым. Также, возможно, рассмотрим другие модификации муравьиного алгоритма, погрузимся чуть глубже, а потом уже перейдем к генетике.
Сейчас предлагаю посмотреть, как безусловный лидер (пока до 100 вершин) идет к глобальному лучшему пути.
Также, предлагаю посмотреть на лидера по времени и количеству вершин SA, который выдает приближенное решение 1000 вершин за 4 000 000 итераций в 34 секунды. Если в прошлой статье 2 000 000 итераций для 500 вершин составляло 90 секунд, то сейчас всего 14.6!
Как-то так. Исходники с комментариями прилагаю. Старался сохранять баланс между читаемостью кода и скоростью. Советую их просмотреть, даже тем кто совсем не знаком с MatLab-ом, так как это сильно поможет вникнуть в суть алгоритмов.
% метод отжига ( на примере задачи коммивояжера )
%--------------------------------------------------------------------------
tic
% clearvars -except cities
clearvars
% -----------------------------ИТЕРАЦИИ------------------------------------
% кол-во итераций
m = 1000000;
% -----------------------------ПАРАМЕТРЫ-----------------------------------
% начальная температура
Tstart = 100000;
% конечная температура
Tend = 0.1;
% начальная температура для вычислений
T = Tstart;
% расстояние
S = inf;
% количество городов
n = 500;
% рисовать графику?
g = 1;
% --------------------------------ПАМЯТЬ-----------------------------------
% матрица расстояний
dist = zeros(n,n);
% -------------------------------------------------------------------------
% генерация городов (x,y)
cities = rand(n,2)*100;
% формируем заранее список случайных чисел
RANDONE = rand(m,1);
% формируем два случайных города заранее
D = randi(n,m,2);
% задаем случайный маршрут
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 = D(k,[1,2]);
% если тут не понятно, посмотрите код из первой части статьи.
if transp(1) < transp(2)
if transp(1) ~= 1 && transp(2) ~= n
S = dist(ROUTE(transp(1)-1),ROUTE(transp(1))) + ...
dist(ROUTE(transp(2)),ROUTE(transp(2)+1));
elseif transp(1) ~= 1 && transp(2) == n
S = dist(ROUTE(transp(1)-1),ROUTE(transp(1))) + ...
dist(ROUTE(transp(2)),ROUTE(1));
elseif transp(1) == 1 && transp(2) ~= n
S = dist(ROUTE(end),ROUTE(transp(1))) + ...
dist(ROUTE(transp(2)),ROUTE(transp(2)+1));
end
else
if transp(2) ~= 1 && transp(1) ~= n
S = dist(ROUTE(transp(2)-1),ROUTE(transp(2))) + ...
dist(ROUTE(transp(1)),ROUTE(transp(1)+1));
elseif transp(2) ~= 1 && transp(1) == n
S = dist(ROUTE(transp(2)-1),ROUTE(transp(2))) + ...
dist(ROUTE(transp(1)),ROUTE(1));
elseif transp(2) == 1 && transp(1) ~= n
S = dist(ROUTE(end),ROUTE(transp(2))) + ...
dist(ROUTE(transp(1)),ROUTE(transp(1)+1));
end
end
%-------------------------------------------------------------------------
if transp(1) < transp(2)
ROUTEp(transp(1):transp(2)) = ROUTEp(transp(2):-1:transp(1));
if transp(1) ~= 1 && transp(2) ~= n
Sp = dist(ROUTEp(transp(1)-1),ROUTEp(transp(1))) + ...
dist(ROUTEp(transp(2)),ROUTEp(transp(2)+1));
elseif transp(1) ~= 1 && transp(2) == n
Sp = dist(ROUTEp(transp(1)-1),ROUTEp(transp(1))) + ...
dist(ROUTEp(transp(2)),ROUTEp(1));
elseif transp(1) == 1 && transp(2) ~= n
Sp = dist(ROUTEp(end),ROUTEp(transp(1))) + ...
dist(ROUTEp(transp(2)),ROUTEp(transp(2)+1));
end
else
ROUTEp(transp(2):transp(1)) = ROUTEp(transp(1):-1:transp(2));
if transp(2) ~= 1 && transp(1) ~= n
Sp = dist(ROUTEp(transp(2)-1),ROUTEp(transp(2))) + ...
dist(ROUTEp(transp(1)),ROUTEp(transp(1)+1));
elseif transp(2) ~= 1 && transp(1) == n
Sp = dist(ROUTEp(transp(2)-1),ROUTEp(transp(2))) + ...
dist(ROUTEp(transp(1)),ROUTEp(1));
elseif transp(2) == 1 && transp(1) ~= n
Sp = dist(ROUTEp(end),ROUTEp(transp(2))) + ...
dist(ROUTEp(transp(1)),ROUTEp(transp(1)+1));
end
end
%--------------------------------------------------------------------------
if Sp < S
ROUTE = ROUTEp;
iter = k;
else
% вычисляем вероятность перехода
P = exp((-(Sp - S)) / T);
if RANDONE(k) <= P
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.')
msgbox ('Выполнено!')
% очищаем переменые
clearvars -except cities ROUTE S iter
% смотрим время
toc
% муравьиный алгоритм ( на примере задачи коммивояжера )
% -------------------------------------------------------------------------
tic
% clearvars -except cities
clearvars
% ------------------------------ИТЕРАЦИИ-----------------------------------
% кол-во итераций ( поколений )
age = 2000;
% кол-во муравьев в поколении
countage = 10;
% кол-во городов
n = 50;
% ------------------------------ПАРМЕТРЫ-----------------------------------
% альфа - коэффициент запаха, при 0 будем ориентироваться только на
% кратчайший путь
a = 1;
% бета - коэффициент расстояния, при 0 будем
% ориентироваться только на оставляемый запах
b = 2;
% коэффициент обновления, глобальное
e = 0.1;
% коэффициент обновления, локальное
p = 0.1;
% количество выпускаемых феромонов
Q = 1;
% баланс между лучшим городом и как в AS
q = 0.9;
% начальный феромон
ph = Q/(n*2000);
% -------------------------------ПАМЯТЬ------------------------------------
% матрица расстояний
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);
% максимальное значение вероятности
val = zeros(1);
% присваем номер города
getcity = zeros(1);
% индекс максимального значения вероятности
indexP = zeros(1);
% максимальное
minDISTiterration = zeros(1);
% -------------------------------------------------------------------------
% генерация городов (x,y)
cities = rand(n,2)*100;
% матрица начальных феромонов
tao = ph*(ones(n,n));
tao(logical(eye(size(tao)))) = 0;
% создаем матрицу расстояний и матрицу обратных расстояний
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;
% смотрим в какой город осуществляется переход
RANDONE = rand;
if RANDONE <= q
[val, getcity] = max(P);
else
% получаем вероятности перехода ( сумма строк должна быть = 1 )
P = P ./ sum(P);
getcity = find(cumsum(P) >= RANDONE, 1, 'first');
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-ых муравьев age-ого поколения
DISTant(k) = S;
% присваевыем лучший маршрут и S
if DISTant(k) < bestDIST
bestDIST = DISTant(k);
bestROUTE = ROUTEant(k,[1:end,1]);
iter = iterration;
end
% вектор "последних" городов k-ых муравьев ( выбирается для старта
% муравьев нового поколения с тех городов, где закончили путь
% предыдущее поколение)
% lastROUTEant = ROUTEant(1:end,end);
% локальное обновление феромона, после каждого муравья
for tL = 1:n
xL = ROUTE(tL);
yL = ROUTE(tL+1);
% считаем новый феромон
tao(xL,yL) = (1-p)*tao(xL,yL) + p*ph;
tao(yL,xL) = (1-p)*tao(yL,xL) + p*ph;
end
end
% --------------------------ГЛОБАЛЬНОЕ ОБНОВЛЕНИЕ--------------------------
% Испаряем феромоны "старого" пути е - коэффициент испарения
tao(tao < 2.500000000000000e-150) = 2.500000000000000e-150;
% для каждого города
for t = 1:n
xG = bestROUTE(t);
yG = bestROUTE(t+1);
% считаем новый феромон
tao(xG,yG) = tao(xG,yG) + e*(Q/bestDIST);
tao(yG,xG) = tao(yG,xG) + e*(Q/bestDIST);
end
end
% строим графику
citiesOP(:,[1,2]) = cities(bestROUTE(:),[1,2]);
plot([citiesOP(:,1);citiesOP(1,1)],[citiesOP(:,2);citiesOP(1,2)],'.r-')
disp (num2str(bestDIST))
msgbox ('Выполнено!')
clearvars -except cities bestDIST bestROUTE iter
toc
Спасибо за внимание. До новых встреч.
Статьи зарубежных авторов:
[1] — M. Dorigo, L.M. Gambardella, Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem // IEEE Transactions on Evolutionary Computation Vol. 1, 1, стр. 53–66, 1997 г.
[2] T. Stützle, H. Hoos, «MAX-MIN Ant System and local search for the traveling salesman problem» // IEEE International Conference on Evolutionary Computation, стр. 309–314, 1997 г.
[3] T. Stützle, M. López-Ibáñez, P. Pellegrini, M. Maur, M. de Oca, M. Birattari, Michael Maur, M. Dorigo, «Parameter Adaptation in Ant Colony Optimization» // Technical Report, IRIDIA, Université Libre de Bruxelles, 2010 г.
