Оптимізація на прикладі. Мурашиний алгоритм (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, для більшого ефекту.

image

Дивимося на метод відпалу при 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.0044939 x^2 + 0.72097 x + 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;
% -----------------------------ПАРАМЕТРИ-----------------------------------
% початкова температура
Сторінку = 100000;

% кінцева температура
Tend = 0.1;

% початкова температура для обчислень
T = Сторінку;

% відстань
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 = Сторінку / 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.500000000000000 e-150) = 2.500000000000000 e-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 р.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.