Метод имитации отжига

    Дорогие друзья, доброго времени суток!

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

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

    Для начала тем, кто забыл, напомню, что же за вид задач мы пытаемся решать такими экзотическими методами. Существует целый ряд задач, чьё точное решение крайне трудно найти в силу большого количества изменяемых параметров. В общем случае у нас есть функция F, чей минимум (максимум) нам нужно найти, и набор её параметров x1..xn, каждый из которых может независимо меняться в своих пределах. Параметры и, соответственно, функция могут изменяться как дискретно, так и непрерывно. Ясно, что у такой функции, скорее всего, будет сложная зависимость от входных данных, имеющая множество локальных минимумов, тогда как мы ищем минимум глобальный. Конечно, всегда остаётся вариант решить задачу полным перебором, но он работает только в дискретном случае и при крайне ограниченном наборе входов. Вот тут и начинают действовать оптимизационные алгоритмы.

    Чтобы не быть сухим, сразу расскажу о задаче, которую я буду решать в рамках этой статьи. Представим себе клетчатую доску размером m на l. Сторона клетки равна единице. На доске нужно найти такую расстановку n точек, чтобы длина минимальной ломаной, соединяющей эти точки, была максимальной. Обычная минимаксная задача. Очевидно, что в данном случае функция F считает длину минимальной ломаной для определённой расстановки точек, а параметрами её являются x и y координаты точек, т.е. 2*n параметров.

    Собственно, алгоритм

    Итак, есть функция, характеризующая систему, и множество координат, на котором эта функция задана. Прежде всего нужно задать начальное состояние системы. Для этого берётся просто любое случайное состояние. Далее на каждом k-ом шаге мы
    1. сравниваем текущее значение F с наилучшим найденным; если текущее значение лучше — меняем глобальное наилучшее
    2. случайным образом генерируем новое состояние; распределение вероятности для него должно зависеть от текущего состояния и текущей температуры
    3. вычисляем значение функции для сгенерированной точки
    4. принимаем или не принимаем сгенерированное состояние в качестве текущего; вероятность этого решения должна зависеть от разности функций сгенерированного и текущего состояний и, конечно, от температуры (чем выше температура, тем больше вероятность принять состояние хуже текущего)
    5. если новое состояние не принято, генерируем другое и повторяем действия с третьего пункта, если принято — переходим к следующей итерации, понизив температуру (но чаще переход к следующему шагу производят в любом случае, чтобы избежать долгого зацикливания)
    Процесс останавливается по достижении определённой температуры. Вещество остыло в точке с минимальной энергией.
    Зачем же нужна такая сложная схема с вероятностями переходов из точки в точку? Почему нельзя просто передвигаться строго от большей энергии к меньшей? Всё дело в уже упоминавшихся локальных минимумах, в которых решение может просто застрять. Чтобы выбраться из них и найти минимум глобальный, необходимо время от времени повышать энергию системы. При этом общая тенденция к поиску наименьшей энергии сохраняется. В этом и состоит суть метода имитации отжига.

    И снова о точках на доске

    Теперь, когда общая схема алгоритма ясна, вернёмся к нашим баранам. Будем реализовывать задачу на java. Для начала опишем доску с точками на ней.
    Доска:
    public class Board {
        private Point[] points;
        private int N;
        private int width;
        private int height;
    
        public Board(int n, int w, int h) {
            points = new Point[n];
            N = n;
            width = w;
            height = h;
        }
    
        public Board(Board board) {
            points = new Point[board.points.length];
            for (int i = 0; i < points.length; i++)
                points[i] = new Point(board.points[i].getX(), board.points[i].getY());
            this.width = board.width;
            this.height = board.height;
        }
        public int getPointsCount() {
            return points.length;
        }
        public int getWidth() {
            return width;
        }
        public int getHeight() {
            return height;
        }
        public Point getPoint(int index) {
            return points[index];
        }
        public void setPoint(int index, Point p) {
            points[index] = p;
        }
    }
    И точка:
    
    public class Point {
        private int x;
        private int y;
    
        public Point(int x, int y) {
            this.x = x;
            this.y = y;
        }
        public double dist(Point p) {
            return Math.sqrt((x - p.x) * (x - p.x) + (y - p.y) * (y - p.y));
        }
        public void move(int dx, int dy, int xBorder, int yBorder){
            if(((x + dx) < xBorder) && ((x + dx) >= 0))
                x += dx;
            if(((y + dy) < yBorder) && ((y + dy) >= 0))
                y += dy;    
        }
        public int getX(){
            return x;
        }
        public int getY(){
            return y;
        }
    }
    Так как для вычисления искомой функции нам важен порядок следования точек, введём их упорядоченный набор:
    
    public class Polyline {
        public Point p[];
        public int N = 5;
    
        public Polyline(int Num) {
            N = Num;
            p = new Point[N];
        }
        public double dist() {
            double d = 0;
            for (int i = 0; i < (N - 1); i++) {
                d += p[i].dist(p[i + 1]);
            }        
            return d;
        }
    }
    Теперь нам нужна функция F, ищущая минимальную ломаную заданной расстановки и считающая её длину.
    public class MinPolyline {
    private double bestDist;
    private Polyline bestMinPolyl;
    private Board board;
    
    public Polyline F(Board b) {
            N = b.getPointsCount();
            board = b;
            int[] perm = new int[N];
            boolean used[] = new boolean[N];
    
            for (int i = 0; i < N; i++) {
                perm[i] = i + 1;
                used[i] = false;
            }
            bestDist = Double.MAX_VALUE;
            permute (0, perm, used);
            return bestMinPolyl;
        }
    
    void permute(int first_index, int[] perm, boolean[] used) {
            if (first_index == N) {
                Polyline polyl = new Polyline(N);
    
                for (int i = 0; i < N; i++){
                    polyl.p[i] = board.getPoint(perm[i] - 1);
                }
                if (bestDist > polyl.dist()){
                    bestDist = polyl.dist();
                    bestMinPolyl = polyl;
                }
                return;
            }
    
            for (int i = 1; i < N+1; i++) {
                if (used[i - 1]) continue;
                perm[first_index] = i;
                used[i - 1] = true;
                permute(first_index + 1, perm, used);
                used[i - 1] = false;
            }
    }
    }

    Ну что ж, все приготовления выполнены. Теперь мы можем начинать реализовывать сам алгоритм.

    Реализация алгоритма

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

    или в коде:
    h = 1 / (1 + Math.exp(-(minPolyl.F(board).dist()- maxDist) / T));
    Таким образом, вероятность принятия позиции с худшим результатом функции (с большей энергией) тем меньше, чем ниже температура и чем больше разница между текущей энергией и оптимальной. Далее установим закон убывания температуры. Отличительной особенностью сверхбыстрого отжига является то, что изменения по всем координатам рассматриваются независимо, и даже температура для каждой координаты определяется собственная:

    Где D — количество координат (т.е. 2*n), k — номер хода. Однако для простоты сделаем всё-таки общую температуру:
    T = To * Math.exp(-decrement * Math.pow(i, (double)(1) / (2*(double)N)));
    Как видно, температура будет убывать экспоненциально, и процедура закончится быстро.
    Наконец, распределение для новой координаты. Следующая величина характеризует относительное изменение одной координаты:

    Где альфа — равномерно распределённая на отрезке [0, 1] величина. Если новая координата не укладывается в рамки своего изменения (в нашем случае — выходит за пределы доски), то расчёт производится снова.
    z = (Math.pow((1 + 1/T), (2 * alpha - 1)) - 1) * T;
    x = board.getPoint(moveNum).getX() + (int)(z * border);
    Теперь у нас есть все необходимые распределения и функции. Осталось лишь собрать всё воедино. В качестве входных данных нужно передать размеры доски, количество точек, начальную и конечную температуру и декремент затухания температурной экспоненты.
    public class AnnealingOptimizer {
    
        public int N;
        public int width;
        public int height;
    
        public AnnealingOptimizer(int n, int width, int height) {
            N = n;
            this.width = width;
            this.height = height;
        }
    
        public Board optimize(double To, double Tend, double decrement){
    
            Board board = new Board(N, width, height);
            Board bestBoard = new Board(N, width, height);
            Random r = new Random();
            double maxDist = 0;
            double T = To;
            double h;
            MinPolyline minPolyl = new MinPolyline();
    
            for (int j = 0; j < N; ++j){
                board.setPoint(j,  new Point(r.nextInt(height), r.nextInt(width)));
                bestBoard.setPoint(j, board.getPoint(j));
            }
    
            int i = 1;
            while (T > Tend){
                for (int moveNum = 0; moveNum < board.getPointsCount(); ++moveNum){
                    NewX(board, moveNum, T, width, true);
                    NewX(board, moveNum, T, height, false);
                }
    
                if (minPolyl.F(board).dist() > maxDist){
                    bestBoard = new Board(board);
                }else{
                    h = 1 / (1 + Math.exp(-(minPolyl.F(board).dist()- maxDist) / T));
                    if (r.nextDouble() > h){
                        board = new Board(bestBoard);
                    }else{
                        bestBoard = new Board(board);
                    }
                }        
    
                maxDist = minPolyl.F(board).dist();
    
                T = To * Math.exp(-decrement * Math.pow(i, (double)(1) / (2*(double)N)));
    
                ++i;
            }
    
            return bestBoard;
        }
    
        private void NewX (Board board, int moveNum, double T, int border, boolean xCoord) {
    
            Random r = new Random();
            int x;
            double z;
            double alpha = r.nextDouble();
    
            z = (Math.pow((1 + 1/T), (2 * alpha - 1)) - 1) * T;
    
            if (xCoord){
                x = board.getPoint(moveNum).getX() + (int)(z * border);
    
                if ((x < border) && (x >= 0)) {
                    board.getPoint(moveNum).move((int)(z * border), 0, width, height);
                } else {
                    NewX(board, moveNum, T, border, xCoord);
                }
            } else {
                x = board.getPoint(moveNum).getY() + (int)(z * border);
    
                if ((x < border) && (x >= 0)) {
                    board.getPoint(moveNum).move(0, (int)(z * border), width, height);
                } else {
                    NewX(board, moveNum, T, border, xCoord);
                }
            }
        }
    }

    Нужно отметить, что некоторые модификации метода имитации отжига дают статистическую гарантию нахождения оптимального решения. Это значит, что при верном выборе параметров алгоритм найдёт оптимальное решение за разумное время без перебора всех входов. Данная реализация находит ответ для доски 12x12 и 5 точек примерно за 15000 итераций. Напомню, что полное число вариантов расстановки равно 12^10. Очевидно, что трудности настройки параметров стоят такого выигрыша. Кстати, метод отжига по быстроте и точности по крайней мере не проигрывает генетическому алгоритму, а чаще всего опережает его.
    Напоследок наглядный результат действия алгоритма.

    Использованная литература:
    А. Лопатин «Метод отжига»
    Поделиться публикацией
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама
    Комментарии 31
    • НЛО прилетело и опубликовало эту надпись здесь
      • +1
        Самый ожидаемый комментарий в топике.
        • +5
          Да имитирует он.
        • +3
          Есть примеры отжига и в биологии.
          Это ренатурация ДНК.
            • +4
              На практике этот метод используется в 100 раз чаще генетических алгоритмов. По крайней мере на топкодеровских марафонах еще никто генетикой не выиграл, а отжиг там весьма распространен.
              • 0
                В проекте я сравнивал работу метода отжига и генетического алгоритма, и генетика даже близко не показывала результатов, сравнимых с отжигом. Хотя не исключаю, что я криво его реализовал.
                • +4
                  Недостаток simulated annealing'а в том, что в холодном состоянии он становится таким же жадным, как скажем градиентный спуск. Это значит, что он за горячую фазу должен попасть в ту область, где локальный максимум будет решением задачи.
                  Генетика — не серебряная пуля, она безусловно медленнее, но она не жадная, что позволяет (и требует) использовать ее для задач с большим количеством локальных максимумов в пространстве решений.
                  За статью большое спасибо.
                  • 0
                    В этом и заключается трудность настройки параметров отжига. Кроме того можно выбрать метод тушения, когда температура убывает линейно.
              • +4
                Цветовая гамма на скрине не очень понравилась сонным глазам =\
                • +1
                  Сонными глазами такой материал лучше не смотреть.

                  Вам теорию объясни, практику распиши, пример приведи, литературу покажи, а вы еще и на иллюстрации ругаетесь!
                  Вы же все мотивацию у авторов убиваете!
                  Или это показатель деградации лени нашей молодежи? Увы и ах.

                  Не понравилась картинка, предложите свою, может автор заменит, если будет согласен с вашей критикой.
                  • +3
                    Вот, банально конечно, но ваши сонные глазки может и не обидятся :)

                  • +1
                    Это же классика: красный на синем.
                  • –12
                    Наверное, не я один открыл статью, чтобы найти тут что-то весёлое? :)
                    • +10
                      Заголовок, безусловно, определил успех поста на Хабре.
                      • +2
                        Блин, я реально полез под кат, чтобы почитать как генетическими алгоритмами можно имитировать недетский отжиг! :) И даже сразу придумал применение — поставить на форуме бота, чтобы тот там отжигал… ;)
                        • +4
                          Вы бы указали где-нибудь в начале, что «имитация отжига» = «simulated annealing». Похоже, есть несколько различных переводов этого термина на русский, и до сегодняшнего дня я не был знаком ни с одним.
                        • +5
                          Хорошая тенденция в последнее время.
                          Все больше хороших статей об алгоритмах.
                          Дейкстра, генерация ландшафтов, несколько статей о генетических алгоритмах, отжиг и т.д.
                          Уже пора говорить что хабр торт?
                          • 0
                            > Кстати, метод отжига по быстроте и точности по крайней мере не проигрывает генетическому алгоритму, а чаще всего опережает его.

                            может быть потому, что ГА ищет локальный оптимум? ) Одним из главных входных параметров ГА является дельта — минимальная степень изменения функции цели на каждой итерации. Задача считается решенной, если F(n)/F(n-1) < d.

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

                            Но, собственно, узнать по дельте, нашли ли мы глобальный или локальный оптимум, нельзя.
                            • +1
                              Гм… мне почему-то очень и очень напоминает Метод Монте Карло? Просто звучит очень и очень похоже.
                              • +2
                                Всё верно. Это вариация (или развитие) метода Монте Карло.
                              • 0
                                Интересно было бы посмотреть на реализацию подбора наиболее вездеходного «велосипеда» — аналогично генетическому алгоритму, ссылку на который недавно постили на хабре.

                                Ну и сравнить, что эффективнее окажется.
                                • +1
                                  Мне кажется, что нельзя метод отжига противопоставиставлять генетическому алгоритму. Потому что, мы легко получим отжиг из ГА просто заменив метод селекции. То есть отжиг — это просто одна из возможных оптимизаций эволюционных алгоритмов.
                                  • 0
                                    ГА, отжиг и алгоритм роя принадлежат к классу естественных алгоритмов. Противопоставлять их можно, потому что часто выбирают ГА тогда, когда не нужно просто из-за популярности.
                                  • НЛО прилетело и опубликовало эту надпись здесь
                                    • +1
                                      Рекомендую:
                                      Clever Algorithms. Nature-Inspired Programming Recipes

                                      Книга описывающая множество интересных алгоритмов, включая Simulated Annealing и Genetic-семейство. К алгоритмам прилагается псевдокод и код на Ruby.

                                      PDFку можно скачать прямо с сайта бесплатно.

                                      • 0
                                        Картинки сдохли, можно перезалить? Спасибо.
                                        • 0
                                          Сделал
                                          • 0
                                            Спасибо ещё раз. Ссылка на pdf'ку кстати тоже сдохла, но её не сложно нагуглить.
                                        • 0
                                          А правильно ли я понимаю, что в формуле
                                          image
                                          Дельта Е и Т должны быть примерно одного порядка?

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