Пользователь
0,0
рейтинг
19 января 2014 в 21:17

Разработка → Введение в оптимизацию. Имитация отжига из песочницы

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

Для понимания статьи Вам понадобятся минимальные навыки программирования и владение математикой на уровне 9 класса средней школы. Статья рассчитана на людей не знакомых с методами оптимизации или только делающих первые шаги в этом направлении.

image



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

Об оптимизации


Для того, чтобы понять, что такое оптимизация, стоит сначала узнать, для чего она нужна. Сразу скажу, что с «задачами оптимизации» в повседневной жизни мы сталкиваемся с завидной регулярностью, не отдавая себе в этом отчёта.

Допустим, мы, получив очередную зарплату, решаем подбить семейный бюджет. Нужно распределить деньги так, чтобы в первую очередь удовлетворить самые важные потребности (еда, плата за электричество, интернет и т.д.), а остаток можно спустить на что-нибудь приоритетом пониже. Купить на все деньги дизайнерский телефон и умереть от голода – не самое лучшее решение, не находите?

Или другой вариант, мы решили отправиться в поход и собираем ранец. Одна проблема – вещей очень много, а ранец у нас не слишком большой (да и сами мы имеем ограниченную грузоподъёмность). Мы хотим взять как можно больше наиболее важных для нас предметов. Если мы идём в горы, прихватить с собой страховочный трос – хорошее решение, взять вместо него любимую PlayStation Vita, а остальное место забить печеньками – наверное, не очень.

Так что же делает оптимизация? Она ищет то самое, хорошее решение. Возможно, не самое лучшее из всех, но точно хорошее.
Теперь чуточку формальнее. Если заглянуть в википедию, можно без труда найти вот такое определение:

Оптимальное (от лат. optimus — наилучшее) решение — решение, которое по тем или иным признакам предпочтительнее других. Следовательно, оптимизация – это способ нахождения оптимального решения.

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

Чтобы свести всё к единому знаменателю, снова заглянем в википедию:
Оптимизация (математика) — нахождение экстремума (минимума или максимума) целевой функции в некоторой области.

Теперь, когда мы более-менее определились, с тем, что делает оптимизация, появляется резонный вопрос – «Как?».
Конечно, мы можем особо не заморачиваться и перебрать все-все решения, а затем сравнить их. В случае с ранцем, нам, наверное, даже это удастся проделать. Но теперь давайте представим, что наш ранец – вовсе не ранец, а целый поезд с множеством вагонов, а перевозим мы разные товары в другой город. Задача не изменилась (загрузить как можно больше ценного), а целевой функцией будет суммарная стоимость товаров. Сколько времени уйдёт у нас на полный перебор решений? Думаю, что много.

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

Не буду томить и скажу, что методов оптимизации существует великое множество. Настолько много, что их принято разбивать на классы, все их я перечислять не буду. Скажу лишь, что «имитация отжига» относится к классу стохастических (от греч. στοχαστικός — «умеющий угадывать») или вероятностных методов. А причём здесь вероятность и как она поможет отыскать решение – чуть ниже.

Имитация отжига


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

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

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

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

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

image

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

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

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

Теперь мы готовы перейти к более строгому описанию метода. Введём обозначения:
Пусть S – множество всех состояний (решений) нашей задачи. Например, для задачи о ранце это будет множество всевозможных наборов вещей.
Пусть – состояние на i-м шаге алгоритма. .
– температура на i-м шаге.

Для того, чтобы использовать имитацию отжига, нам понадобится определить три функции:
  • Функцию энергии или, проще говоря, то что мы оптимизируем.

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

  • Функцию, порождающую новое состояние.


Теперь остановимся на каждой подробнее.
E каждому решению по какому-то правилу ставит в соответствие число. Зависит от конкретной задачи.

T ставит номеру итерации i в соответствие температуру. Функция определяет, как долго будет работать наш алгоритм (на самом деле она имеет более глубокое значение, но нам это знать пока не обязательно). Если T будет линейной функцией, то время работы будет относительно большим. В случае же со степенной функцией, скажем , всё закончится очень быстро, но далеко не факт, что мы не застрянем в низменности, так не добравшись до финиша.

Для этой статьи я выбрал самый простой, линейный вариант.

function [ T ] = DecreaseTemperature( initialTemperature, i)
%initialTemperature - начальная температура
%i - номер итерации
T = initialTemperature * 0.1 / i; 
end


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

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

Теперь мы можем в самом общем виде описать алгоритм метода имитации отжига:
  • На входе: минимальная температура , начальная температура
  • Задаём произвольное первое состояние
  • Пока
    • Если , тогда
    • Если переход осуществляется с вероятностью
    • Понижаем температуру

  • Возвращаем последнее состояние s


Как видите, всё очень просто, отдельно стоит пояснить только переход в новое состояние. Если энергия «кандидата» меньше, он становится новым состоянием, в противном случае, переход будет вероятностным (поэтому метод относят к классу стохастических).

Чтобы это точно не вызывало трудностей, я позволю себе напомнить некоторые самые базовые понятия теории вероятностей.

Неформально определим вероятность, как меру возможности наступления некоторого события. Вероятность достоверного события (событие, которое точно случится) равна единице, невозможного – нулю. Все остальные значения находятся между.
Например, при броске игральной кости вероятность выпадения любого числа равна единице, а тройки – 1/6.
В нашем случае, подставив в формулу значения разницы энергий и температуры мы получим вероятность перехода. Остаётся только этот переход осуществить, но как это сделать?

Код вычисления вероятности:
function [ P ] = GetTransitionProbability( dE, T )
     P = exp(-dE/T);
end

Давайте представим себе единичный отрезок и разобьём его на две части относительно значения нашей вероятности.



Затем воспользуемся функцией генерации равномерно распределенного случайного числа, в Matlab это rand (как и почти во всех остальных языках). Теперь отметим это число на отрезке.



Если число попало в левую часть – осуществляем переход, в правую – не делаем ничего.
function [ a ] = IsTransition(probability)
    value = rand(1);

    if(value <= probability)
        a = 1;
    else
        a = 0; 
    end

end


Теперь, думаю, всё понятно, и мы можем перейти непосредственно к решению нашей первой задачи.

Задача коммивояжёра


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

Для начала давайте определимся с данными. Пусть наши города случайным образом разбросаны в квадрате 10x10. Каждый город, соответственно, представлен парой координат. Всего 100 городов.


Данные можно без труда получить в Matlab при помощи следующего кода:
data = rand(2,100)*10;

В других языках придётся написать цикл, но по сути это ничего не меняет. Мы генерируем случайные числа уже знакомой нам функцией rand() и домножаем их на 10, чтобы они лежали в заданной области.

Решение задачи методом имитации отжига


Обозначим множество всех городов C, общее количество обозначается |C| и в нашем случае равно 100. Каждый город представляется как пара координат с соответствующим индексом.

По условию, решением является маршрут между всеми городами, значит множество состояний S – это все возможные маршруты, проходящие через каждый город. Другими словами множество всех упорядоченных последовательностей элементов C, в которой всякий встречается ровно один раз. Очевидно, что длина каждой такой последовательности |C|.

Как мы помним, для того, чтобы использовать метод имитации отжига мы должны определить две функции, зависящие от каждой конкретной задачи. Это функция энергии E (или «целевая функция» в общепринятой терминологии) и функция F, порождающая новое состояние.
Так как мы стремимся минимизировать расстояние, оно и будет «энергией». Следовательно, наша целевая функция будет иметь такой вид:



Здесь написано не что иное, как сумма Евклидовых расстояний между парами городов в маршруте . Так как задача замкнутая, мы были вынуждены подписать в конце, это расстояние между последним и первым городом. Если не углубляться, то Евклидово – это то самое расстояние между точками, которое мы все проходили на уроках геометрии. Задаётся оно следующей формулой:



Где

Теперь давайте подумаем, как же нам получать новое состояние? Первое что приходит в голову – менять местами два произвольных города в маршруте. Идея правильная, но такое изменение непредсказуемо влияет на , метод будет работать долго и не факт, что успешно.

Хорошим вариантом в данном случае будет выбрать два произвольных города в маршруте и инвертировать путь между ними. Например, у нас был маршрут (1,2,3,4,5,6,7). Генератор случайных чисел выбрал города 2 и 7, мы выполнили процедуру и получили (1,7,6,5,4,3,2).

Ниже представлен код для функции F:
function [ seq ] = GenerateStateCandidate(seq)
	%seq - предыдущее состояние (маршрут), из которого 
    %мы хотим получить состояние-кандидат
    
    n = size(seq,1); % определяем размер последовательности
    i = randi(n,1); % генерируем целое случайное число
    j = randi(n, 1);% генерируем целое случайное число
        
    if(i > j) 
        seq(j:i) = flipud(seq(j:i)); % обращаем подпоследовательность
    else
        seq(i:j) = flipud(seq(i:j));% обращаем подпоследовательность
    end

end


Теперь у нас есть всё что нужно, и мы можем начать процесс оптимизации. Но прежде чем мы перейдём к самому интересному, я приведу полный код программы. Думаю, будет весьма занятно сравнить его с приведённым парой абзацев выше псевдокодом:
function [ state ] = SimulatedAnnealing( cities, initialTemperature, endTemperature)

    n = size(cities,1); % получаем размер вектора городов

    state = randperm(n); % задаём начальное состояние, как случайный маршрут
    % Функция randperm(n) - генерирует случайныую последовательность из целых чисел от 1 до n
    
    currentEnergy = CalculateEnergy(state, cities); % вычисляем энергию для первого состояния
    T = initialTemperature;
    
    for i = 1:100000  %на всякий случай ограничеваем количество итераций
    % может быть полезно при тестировании сложных функций изменения температуры T       

        stateCandidate = GenerateStateCandidate(state); % получаем состояние-кандидат
        candidateEnergy = CalculateEnergy(stateCandidate, cities); % вычисляем его энергию
        
        if(candidateEnergy < currentEnergy) % если кандидат обладает меньшей энергией
            currentEnergy = candidateEnergy; % то оно становится текущим состоянием
            state = stateCandidate;
        else
            p = GetTransitionProbability(candidateEnergy-currentEnergy, T); % иначе, считаем вероятность
            if (MakeTransit(p)) % и смотрим, осуществится ли переход
                currentEnergy = candidateEnergy;
                state = stateCandidate;
            end
        end;

        T = DecreaseTemperature(initialTemperature, i) ; % уменьшаем температуру
        
        if(T <= endTemperature) % условие выхода
            break;
        end;
    end

end



Как видите, ничего сложного.

Результаты


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

Программа запускалась с = 10 и = 0.00001.

На рисунке изображено состояние на самой первой итерации. Вряд ли коммивояжер испытал бы огромную радость предложи мы ему такой маршрут.


Результат на выходе выглядит куда приятнее:


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

Исходные коды можно скачать здесь.

Вместо заключения


Понимаю, что было бы интересно взглянуть на сравнение «имитации отжига» с другими методами оптимизации, например, генетическими алгоритмами или динамическим программированием; привести примеры из интересных областей вроде machine learning; указать на связь с Марковскими моделями и много чего ещё. Однако, в рамках одной статьи всего не объять.

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

P.S.

Пользуясь случаем, хотел бы поблагодарить Антонову Анну и Дмитрия Игнатова, которые прочли статью перед публикацией и дали мне несколько хороших советов.
@DukeGonzo
карма
54,2
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама

Самое читаемое Разработка

Комментарии (36)

  • +2
    Сколько итераций было сделано для получения результата на второй картинке?
    • 0
      Для данной функции T выходит ровно 100 000.
  • +2
    Не могли бы вы добавить в статью сравнение этого метода с генетическим алгоритмом? Плюсы и минусы этих методов относительно друг-друга.
    Ведь в целом они похожи, но хотелось бы увидеть, когда какой следует применять.
    Я вот несколько раз применял ген. алг для различных задач, в частности для обучения многослойного перцептрона — получил весьма интересные результаты, но с методом из вашей статьи знаком только понаслышке и не знаю его преимуществ и недостатков относительно ген. алга (в частности, в вопросе обучения НС)

    Upd: не увидел ваше заключение, прошу прощения. Но вопрос остается открытым, и, на мой взгляд, ответ на него даже важнее реализации метода имитации отжига в матлабе.
    • +1
      Самому было бы интересно досконально разораться. Подумываю, написать отдельную статью об этом.
  • +2
    Большое спасибо за интересную статью!
    С радостью почитаю обзор различных методов оптимизации. Сам занимаюсь целочисленным программированием (точнее 0-1), но, боюсь, фундамента моих знаний недостаточно для столь понятного и отличного описания. Оптимизация — очень широкая область, но на Хабре, к сожалению, почти нет материалов по ней.
  • +2
    Интересный метод. Спасибо за наглядное описание.
    Надо бы попробовать его приспособить к задачам решения нелинейных уравнений. Иногда, вместо решения сложной системы многих переменных f(x_1..x_n) = 0 можно попытаться минимизировать каким-нибудь образом построенную из этих уравнений функцию, например, сумму их квадратов (что абсолютно естественно для комплексных чисел — минимизация квадрата модуля). Возможно, этот метод себя сможет хорошо проявить, в отличие от многомерных секущих (плохо сходятся) или Ньютона-Рафсона (жутко долго писать и неудобно вычислять производные, если функция не задана аналитически).
  • +4
    Следует также отметить, что метод имитации отжига — это слега адаптированный алгоритм Метрополиса (Metropolis algorithm), который, в свою очередь, является примером метода Markov Chain Monte Carlo (MCMC). А MCMC, по мнению нескольких научных изданий, входит в 10-ку наиболее важных для инженерии алгоритмов 20-ого столетия (для примера, туда же обычно включают быструю сортировку).

    Спасибо за прекрасное описание! Жирный плюс в карму, и будет ждать следующих статей!
    • +2
      Очень полезное дополнение, особенно с учетом того, что MCMC почти не освещался на Хабре. Всем для первого знакомства с темой рекомендую вот такую книгу, написано очень просто, но в то же время достаточно строго.
      • 0
        Не знаете, эта книжка есть на русском?
  • 0
    Супер!
    Поиск кратчайшего пути по «бразильской» системе!
    Очень напоминает как много-много стартапов 99% из которых проваливается в итоге рождают что-нибудь в духе эппл или макрософт!
    Конечно в реальной жизни множество особенностей, но метод отжига заполнил очень большой кусок в картине того, как реальный человек (или структура состоящая из людей) взаимодействует с миром!
    И прямо так руки к чешутся взять и применить этот алгоритм к какой-нибудь олимпиадной задачке из школьных времен. Такое чувство, что множество хитрых задачек сведутся к одному алгоритму и будут проходить 90% тестов :)
    • +1
      При словосочетании «олимпиадная задача» надолго завис, размышляя о задаче траты средств на олимпиаду 2014 в Сочи. Ох, что за времена.
    • 0
      Такое чувство, что множество хитрых задачек сведутся к одному алгоритму и будут проходить 90% тестов :)


      Только вот скорее всего (в олимпиадах, проводящихся по правилам ACM) за это дадут 0 баллов :) Нужно пройти все тесты.
  • +4
    А о самом интересном-то и не написали. Почему вероятность перехода именно exp(-dE/T)?
    • +10
      Собственно, в этом и сидит физический аспект алгоритма.
      Чистая термодинамика и никакой магии.
      Вероятность пребывания системы в состоянии с энергией E при температуре T равна C*exp(-E/T), где С — нормировочный множитель. Это свойство общее для любых макроскопических систем и известно под названием распределения Гиббса. Вероятность состояния с энергией повыше равна, соответственно, C*exp(-(E+dE)/T).

      А вероятность перехода между двумя этими состояниями равна экспоненте, в которой вместо энергии стоит dE — если условно принять нижний уровень за основное состояние (считать, что для него E=0 — в случае с потенциальной энергией именно это и можно сделать), как раз и получим, что вероятность попасть в состояние с энергией dE равна exp(-dE/T), только нормировка C здесь будет пересчитана.
  • 0
    Очень интересная и крайне доступно написанная статья! Спасибо!
    Что можете сказать про скорость такого алгоритма и точность относительно других методов?
    Я занимался подобными изысканиями для поиска лучших торговых стратегий в трейдинге. Мне нужно было найти не только глобальный экстремум, но и максимально исследовать все локальные экстремумы! Разработал новый алгоритм, который не стремится к глобальному экстремуму, а удаляет из пространства исследования все минимумы и сходится ко всем экстремумам.
    Метод оказался эффективнее и удобнее по сравнению с теми, которые я ранее использовал. Буду рад, если вы как эксперт оцените мой алгоритм, а то я прямых аналогов ему пока не нашел.
    Удачи!
  • +2
    Понимаю, что было бы интересно взглянуть на сравнение «имитации отжига» с другими методами оптимизации, например, генетическими алгоритмами или динамическим программированием; привести примеры из интересных областей вроде machine learning; указать на связь с Марковскими моделями и много чего ещё. Однако, в рамках одной статьи всего не объять.


    жду продолжения
  • +1
    Спасибо, отличная статья!
  • 0
    Спасибо за статью.
    А можно выложить куда-то полные исходники?

  • +1
    Спасибо, статья пришлась кстати.
    Попробую «метод отжига» (MO) при случае на своих реальных данных.

    На практике приходится работать с оптимизацией и еще один метод не помешает.

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

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

    Так, например, может и не быть большого смысла применять МО в задаче коммивояжера с простым (как в примере) — равномерным распределением городов по некоторой плоскости.

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

    Другая возможность — если функция генерации нового состояния почему-то из каких-то тактических соображений работает не затрагивает весь набор «городов.

    Вид и свойствацелевой функции (как вычисляется „энергия“) тоже будут влиять.
    А если без извлечения корня?...
    Вот тут, кстати, возникает вопрос… А не лучше ли будет обойтись без извлечения корня?..
    Сумма квадратов разностей имеет более „энергетический“ вид с точки зрения меня, как физика…
    Ок, еще один вопрос для „поиграться на досуге.


    Захотелось поиграться с МО самому.
    Начал с ММС, чтоб потом сравнить.
    Вот здесь под спойлер положил на скорую руку сооруженный R код для МНС и пару картинок.
    R код для монотонного спуска, для любителей
    Может кого заинтересует, например, в порядке знакомства с R.
    А то я все собираюсь написать пост с примерами практического использования.
    Заодно можно сравнить с авторским матлабовским кодом.

    Код не идеален, но рабочий.
    ####### Init
    n.cities=100
    
    ## Матрица координат
    coors=data.frame(cbind(x=runif(n.cities), y=runif(n.cities)))
    
    ## Случайно выбранный маршрут -- вектор индексов городов в случайной последовательности
    s0=sample(1:n.cities)
    
    ######## Functions ##############
    plot.s<-function(s=s0) {
    	crs=rbind(coors[s,], coors[s[1],])
    	plot(crs$x,crs$y,type="b",axes=F)
    }
    
    E<-function(s=s0){
    	crs=rbind(coors[s,], coors[s[1],])
    	sum(sqrt(diff(crs$x)^2+diff(crs$y)^2))	
    }
    
    	
    next.s<-function(s=s0){
    	ii=sample(1:length(s),2)
    	s[ii[1]:ii[2]]<-s[ii[2]:ii[1]]
    	return(s)
    }
    
    runOnce<-function(s=s0, n.iter=100000, plo=T) {
    	Es=E(s)
    	for (i in 1:n.iter){
    		ns=next.s(s)
    		Ens=E(ns)
    		cat(i, Es,Ens,"\n")
    		if (Ens<Es){
    			s=ns
    			Es=Ens
    			if (plo){
    				plot.s(s)
    				title(paste("i=",i,"E=",signif(Es,3)))
    			}
    		}
    	}
    	return(s)
    }
    
    ######## End of Functions #######
    
    ## Main ##
    plot.s(s0)
    title(paste("Randomly generated route,",n.cities," cities, E=",signif(E(s0),3)))
    
    ### Один "сеанс"
    
    runOnce()
    
    ## Результаты 16 попыток
    op<-par(mfcol=c(4,4),mar=c(0,0,2,0))
    for (j in 1:16){
    	s=runOnce(plo=F,n.iter=100000)
    	plot.s(s)
    	title(paste("E=",signif(E(s),3)))	
    }
    

    Руки не дошли рассмотреть с “отжигом».

    Картинки:
    «Города» и случайный маршрут обхода.
    В шапке — значение «энергии» — полной длины данного маршрута.
    image

    16 попыток
    image


    Не совсем прозрачно пока для меня вопрос «температуры».
    Выдастся время на еще один подзод — попробую довести реализацию своего кода до МО.

    Будем думать.
    • +2
      Почему-то не видны картинки, залитые на hfbrastorage…
      Перевставил.
      Скрытый текст
      Картинки:
      «Города» и случайный маршрут обхода.
      В шапке — значение «энергии» — полной длины данного маршрута.
      image

      16 попыток, монотонный спуск.
      image

      • +1
        Мне по работе нужно к своим данным подобрать подходящий метод оптимизации, поэтому я меркантильно заинтересован, чтоб дискуссия не затухала. Тема интересует многих — новые статьи появляются
        Скрытый текст

        (правда, качеством пониже чем эта, не в обиду другим авторам).


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

        Выложил свой код R-код для задачи коммивояжера.
        Пока без отжига, поскольку еще в процессе экспериментов с поведением температуры.

        Заодно вожусь с ГА, чтобы сравнить с «отжигом». Может позже подискутируем и об этом, если найдутся желающие.

        Нижеследующий код производит расчет и строит картинки (см. мой комментарии выше).

        Один прогон из 100к итераций совершается в R 3.01 около 1.5 мин на маке c i5.

        n.cities=100
        
        ####### Init
        ## Матрица координат
        coors=data.frame(cbind(x=runif(n.cities), y=runif(n.cities)))
        
        ## Случайно выбранный маршрут -- вектор индексов городов в случайной последовательности
        s0=sample(1:n.cities)
        
        ######## Functions ##############
        
        # Нарисовать один вариант
        plot.s<-function(s=s0) {
        	crs=rbind(coors[s,], coors[s[1],])
        	plot(crs$x,crs$y,type="b",axes=F)
        }
        
        #Подсчитать длину кольцевого пути
        E<-function(s=s0){
        	crs=rbind(coors[s,], coors[s[1],])
        	sum(sqrt(diff(crs$x)^2+diff(crs$y)^2))	
        }
        
        
        #Сгенерировать следующий вариант рутем изменения направления обхода между двумя выбранными городами на противоположный	
        next.s<-function(s=s0){
        	ii=sample(1:length(s),2)
        	s[ii[1]:ii[2]]<-s[ii[2]:ii[1]]
        	return(s)
        }
        
        # Зарустить одн раз
        runOnce<-function(s=s0, n.iter=100000, plo=T) {
        	Es=E(s)
        	for (i in 1:n.iter){
        		ns=next.s(s)
        		Ens=E(ns)
        		cat(i, Es,Ens,"\n")
        		if (Ens<Es){
        			s=ns
        			Es=Ens
        			if (plo){
        				plot.s(s)
        				title(paste("i=",i,"E=",signif(Es,3)))
        			}
        		}
        	}
        	return(s)
        }
        
        ######## End of Functions #######
        
        #########################
        ## Main 
        ## Здесь лучше работать в интерактивном режиме
        plot.s(s0)
        title(paste("Randomly generated route,",n.cities," cities, E=",signif(E(s0),3)))
        
        ### Один "сеанс"
        st=system.time(runOnce())
        mtext(paste("  time(s): ",st["elapsed"]))
        
        ## Результаты 16 попыток
        op<-par(mfcol=c(4,4),mar=c(0,0,2,0))
        for (j in 1:16){
        	s=runOnce(plo=F,n.iter=100000)
        	plot.s(s)
        	title(paste("E=",signif(E(s),3)))	
        }
        
        


  • 0
    Вообще почти ничего не смыслю в данной области, поэтому возможно скажу какую-то очевидную вещь. Но мне кажется алгоритм наверное можно было бы ускорить, если ввести некую оценочную фукцию, которая бы отражала разницу энергий между итерациями. Как только оценочная фукция перестает уменьшаться (или перестает достаточно быстро или стабильно уменьшаться) можно прекращать итерации.
    • 0
      Так это локальный минимум может быть. Достаточно широкое «плато».
      • 0
        ну да, полочка должна быть довольна широкая. Ну ее ширину у локальных минимумов можно определить эмперически.
  • 0
    Спасибо за статью, в свое время, когда учился, использовал для решения задачи коммивояжёра генетические алгоритмы.
    Из замечаний, в статье отсутствует реализация функции CalculateEnergy().
    Не могли бы вы ее дополнить?
  • +1
    Проясните пожалуйста момент с генерацией нового расстояния. Почему "менять местами два произвольных города в маршруте" меньше влияет на E, чем "инвертировать путь".

    В первом случае мы видим изменение расстояний в четырёх местах
    1*7*3,4,5,6*2*

    А во втором, только в двух (ведь если мы идём в обратную сторону, расстояние не меняется)
    1*7,6,5,4,3,2*
  • +3
    По личному опыту, на практике этот алгоритм работает лучше генетического. Почему так с точки зрения теории — сказать не возьмусь. Люди говорят, что в целом на практике отжиг применяется гораздо чаще.

    Думаю, надо отметить более акцентировано, что трудность использования метода имитации отжига заключается именно в настройке функции температуры.

    P.S. Я тоже писал об этом алгоритме.
    • 0
      Да, кстати, спасибо и за ту статью.
      Меня немного отпугнула Java, поскольку я уже больше 20 лет, как сам отошел к скриптовым языкам.

      Сейчас посмотрел ее — кое что для себя «выцепил» новое.
      Очень полезными нашел комментарии к той статье.

      В своих задачах я как-то применял ГА и меня устраивало.

      Сейчас у меня назревает следующий раунд в связи с поступлением новых данных. И я, похоже, запАл на метод отжига.
      Так что ждите реализации в R.
  • +2
    >Первое что приходит в голову – менять местами два произвольных города в маршруте. Идея правильная, но такое изменение слишком слабо влияет на E, метод будет работать долго и не факт, что успешно.

    По-моему, наоборот. Такое изменение, напротив, сильно влияет на E и, вероятнее всего, его увеличивает.
    Вот что говорит на этот счет английская википедия:

    >In the traveling salesman problem above, for example, swapping two consecutive cities in a low-energy tour is expected to have a modest effect on its energy (length); whereas swapping two arbitrary cities is far more likely to increase its length than to decrease it. Thus, the consecutive-swap neighbour generator is expected to perform better than the arbitrary-swap one
    • 0
      Да, для меня такой способ генерации нового состояния тоже не показался совсем очевидным.
      От свойств данных будет очень сильно все зависеть.

      Но всему критерий практика.
      Иногда быстрее реализовать разные варианты алгоритма и попробовать на своих (или тестовых) наборах данных.

      Я так обычно делаю.
      • +2
        Я не про «неочевидность» говорю, с этим как раз вопросов нет, вполне логично, что новые изменения не должны вызывать совершенно случайное изменение энергии, особенно с понижением температуры, потому что тогда будет под вопросом сходимость алгоритма.

        Я говорю про формулировку в статье, написано, что «такое изменение слишком слабо влияет на энергию», в то время как на практике — совсем наоборот, поэтому его и не используют.
    • +1
      Да, Вы правы, тут ошибка.

      Попробую объяснить почему.
      Если представить города в виде графа, то видно, что для того, чтобы поменять местами два произвольных города нужно изменить четыре ребра. В случае, когда мы «переворачиваем» подпоследовательность — изменяется только два. Следовательно, влияние на E меньше (что хорошо).
      Тогда встаёт вопрос почему вторая эвристика лучше?
      Тут Википедия совершенно права, при перестановке двух произвольных элементов вероятность получить плохое состояние выше. Попытаюсь проиллюстрировать на примере:

      Меняем A и B. Даже получив выигрыш за счёт ребер XB и AY, мы много потеряли на последовательности между A и B (рёбра BZ и WA).


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


      Такое объяснение, конечно, не подходит на роль полноценного доказательства. Постараюсь, когда будет больше свободного времени найти правильные слова и дополнить статью.
      Спасибо за замечание.

  • 0
    В/К статье есть еще несколько вопросов, которые не грех было бы обсудить.

    Трудновато это делать в формате комментария, а на статью-отклик не хватает, к сожалению, прочих ресурсов.
  • 0
    1.Мы случайно экспериментируем с функцией T (временем или количеством итераций, что у вас там). А вдруг при 10T найдется лучшее решение? И при 100T еще лучшее? Как оценить точность приближения? (тест в студию, где заранее полностью просчитаны глобальный минимум и локальные минимумы и при различных T как успешно ваша программа находит глобальный минимум)
    2.Область обсчета и точки заведомо конечны. Что будете делать, когда они огромны (почти бесконечны)? Как оценивать будете, что в таком-то направлении считать невыгодно, а в другом выгодно и в какой момент программа остановится?
    3.Предлагаю вам написать ваш программу на Lua и посмотреть, насколько выгодно она находит глобальный минимум в самых жестких условиях, где ваш алгоритм будет конкурировать с другими известными алгоритмами при поиске максимально выгодного свертывания белков — на foldit.
  • +1
    Очень хорошая и интересная статья!
    Местами, на мой вкус, слишком формализовано, но это я говорю как человек который более-менее «в теме». Для людей, впервые знакомящихся с данным материалом, зато будет легко вникнуть и разобраться.
    Данный пост сподвиг меня на написание своей статьи на эту же тему. Вернее, конкретно эта фраза: «Первое что приходит в голову – менять местами два произвольных города в маршруте. Идея правильная, но такое изменение непредсказуемо влияет на, метод будет работать долго и не факт, что успешно.». Я как раз недавно сам решал задачу коммивояжера и придумал очень похожий на метод отжига алгоритм и использовал как раз идею случайного обмена местами двух городов. Кому интересно — почитайте.

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.