Как я изобретал метод имитации отжига

Доброго времени, Хабр!

Сподвигло меня на написание этой работы прочтение «Введение в оптимизацию. Имитация отжига». Так уж сложилось, что как раз недавно я столкнулся с задачей коммивояжера и для ее решения придумал алгоритм, суть которого, как оказалось, очень близка к идее алгоритма имитации отжига, описываемого в указанной статье. Более того, там даже есть «отсылка» к моей идее, а еще похожие обсуждения велись в комментариях, потому я решил, что сообществу будет интересно посмотреть на реализацию.

Начну с небольшой вводной. Я студент 1-го курса магистратуры на специальности «Программная инженерия», и в этом семестре у нас был курс распределенных систем, в котором мне досталась на реализацию задача коммивояжера. Причем ладно бы так, но еще и с требованием реализовать ее параллельно и с использованием четырех разных технологий (Windows Azure, библиотека openmpi, язык программирования Go, язык программирования Limbo в OS Inferno). Натолкнись я тогда на настолько качественно описанный вариант решения, проблем было бы гораздо меньше, но на тот момент я сумел обнаружить только метод ветвей и границ. Мерзкая штука надо сказать, потому идея реализовывать его параллельно, да к тому же с использованием неизвестных ранее технологий сразу показалась мне неудачной. В результате долгих обдумываний данной задачи я пришел к алгоритму, который оказался на удивление простым, хотя и не очень эффективным. Тем не менее, я его сумел сдать преподавателю, а проведенная «исследовательская» работа над алгоритмом, как оказалось, вполне сгодится для данной статьи. Разве что для статьи я приведу программный код на чистом C# и без использования параллельных заморочек.
Подробно на вопросах что такое оптимизация и метод отжига я останавливаться не буду, считаю, что в статье «Введение в оптимизацию. Имитация отжига» оно объяснено просто замечательно, и всем желающим разобраться подробнее я советую ее прочесть. Я же перейду сразу к описанию задачи коммивояжера и моему решению.

Анализ задачи

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

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

Алгоритм решения

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

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

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

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

Программная реализация на C#

class CCities
    {
        //массив коррдинат городов
        public Point[] Coordinate;
        
        public CCities(int N, int maxValue) //maxValue - размер элемента pictureBox на форме
        {
            Random random = new Random();

            Coordinate = new Point[N];

            //создаем более узкие границы, чем сам pictureBox, чтобы города не лежали с самого краю
            //так просто визуально приятнее выглядит
            int minBorder = (int)(maxValue * 0.05);
            int maxBorder = (int)(maxValue * 0.95);

            for (int i = 0; i < N; i++)
            {
                Coordinate[i] = new Point(random.Next(minBorder, maxBorder),
                                          random.Next(minBorder, maxBorder));
            }
        }
    }

    class CPath
    {
        //расстояния между городами
        double[,] distance;
        //индексы городов формируют искомый путь
        public int[] Path;

        public CPath(CCities map)
        {
            //на вход передаем уже созданные города
            distance = new double[map.Coordinate.Length, map.Coordinate.Length];
            //формируем матрицц расстояний, работать в дальнейшем будем именно с ней
            for (int j = 0; j < map.Coordinate.Length; j++)
            {
                distance[j, j] = 0;

                for (int i = 0; i < map.Coordinate.Length; i++)
                {
                    double value = Math.Sqrt(Math.Pow(map.Coordinate[i].X - map.Coordinate[j].X, 2) +
                                             Math.Pow(map.Coordinate[i].Y - map.Coordinate[j].Y, 2));
                    distance[i, j] = distance[j, i] = value;
                }
            }

            //создаем начальный путь
            //массив на 1 больше кол-ва городов, а первый и последний индексы равны 0 - это сделано для того чтобы "замкнуть" путь
            Path = new int[map.Coordinate.Length + 1];
            for (int i = 0; i < map.Coordinate.Length; i++)
            {
                Path[i] = i;
            }
            Path[map.Coordinate.Length] = 0;
        }

        //метод, реулизующий алгоритм поиска оптимального пути
        public void FindBestPath()
        {
            Random random = new Random();

            for (int fails = 0, F = Path.Length * Path.Length; fails < F; )
            {
                //выбираем два случайных города
                //первый и последний индексы не трогаем
                int p1 = 0, p2 = 0;
                while (p1 == p2)
                {
                    p1 = random.Next(1, Path.Length - 1);
                    p2 = random.Next(1, Path.Length - 1);
                }
                
                //проверка расстояний
                double sum1 = distance[Path[p1 - 1], Path[p1]] + distance[Path[p1], Path[p1 + 1]] +
                              distance[Path[p2 - 1], Path[p2]] + distance[Path[p2], Path[p2 + 1]];
                double sum2 = distance[Path[p1 - 1], Path[p2]] + distance[Path[p2], Path[p1 + 1]] +
                              distance[Path[p2 - 1], Path[p1]] + distance[Path[p1], Path[p2 + 1]];

                if (sum2 < sum1)
                {
                    int temp = Path[p1];
                    Path[p1] = Path[p2];
                    Path[p2] = temp;
                }
                else
                {
                    fails++;
                }
            }
        }

        //возвращает длину пути
        public double PathLength()
        {
            double pathSum = 0;
            for (int i = 0; i < Path.Length - 1; i++)
            {
                pathSum += distance[Path[i], Path[i + 1]];
            }
            return pathSum;
        }
    }

Примеры работы алгоритма

Анализ работоспособности алгоритма

Докажу вначале эффективность распараллеливания такого алгоритма. Так как каждый из параллельных процессов начинает работу с одинаковыми начальными условиями, вероятность найти оптимальное решение для каждого из них будет одинаковой. Примем ее за p. Тогда вероятность не найти такое решение будет равна q = 1 – p. Таким образом, вероятность P того, что хотя бы один из процессов найдет оптимальный путь равна:

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

Можно видеть, что уже при K = 6..7 вероятность P становится равна 0,5. То есть польза от использования нескольких процессов очевидна – очень быстрое увеличение вероятности.
Найдем теперь вероятность p для одного процесса. Очевидно, она обратно пропорциональна размеру N матрицы и прямо пропорциональна значению F, которого достигает счетчик, когда алгоритм прекращает свою работу. Более точную оценку, к сожалению, дать практически невозможно. Чуть ли не единственное, что мы можем сказать о графе – это что количество разных путей равно (N-1)!⁄2 (вычитаем единицу потому что один город является стартовым и в переборе явно не участвует, а делим на два из-за того, что один и тот же путь можно обойти в обе стороны).
Таким образом, оценить искомую вероятность можно лишь экспериментальным путем. Была проведена серия опытов для разных параметров N. В качестве ограничения для счетчика принято значение F = N^2. Для каждого N было проведено 10000 прогонок и полученные результаты занесены в таблицу:

Поясню введенные обозначения:
P(min) – вероятность найти оптимальное решение. P(Х%) – вероятность найти решение, отличающее от оптимального не более чем на Х%. То есть, например, для P(10%) такими решениями будут те, которые попадают в интервал [min, min+(max-min)*0.1], где min и max – минимальный и максимальный соответственно найденные пути.
Тут следует сделать лирическое отступление и сказать вот что. Разумеется, такой анализ не претендует на сколько-нибудь действительную точность, т. к. даже 10000 экспериментов очень мало для получения искомых закономерностей (при повторном запуске результат может оказаться совершенно другим), да и, самое главное, нет гарантии что алгоритм нашел и правда оптимальный из возможных путей. Для оценки он оперирует знанием о лучшем из найденных, но далеко не факт, что он к тому же лучший из всех. Тем не менее, я считаю можно пользоваться полученными результатами пусть не для количественного, но хотя бы для качественного понимания работы алгоритма.
Проиллюстрируем эти же результаты:

Достоинства метода

• «Идейная» простота алгоритма;
• Простота реализации как в последовательном, так и в параллельном вариантах;
• Быстро растущая эффективность алгоритма при его параллельной реализации;
• Высокая вероятность нахождения решения, близкого к оптимальному.

Недостатки метода

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

Оптимизация алгоритма
В процессе работы и исследования метода, я предположил следующие возможные пути улучшения алгоритма:

• Добавление «истории» уже проверенных путей. На данный момент алгоритм не запоминает проверенные пути и никак не контролирует получаемые на каждой итерации новые варианты. Добавление функции проверки и отсеивания таких путей очень сильно должен сказаться на точности результата;
• Возможно, есть способ более «тонкой» настройки контрольного значения счетчика. Как видно из проверочных графиков, при текущих настройках (F = N^2) появляется «яма» в районе матриц размера N = 20. У меня нет точного понимания причины такого поведения графика, но, вероятно, есть способ подбирать это значение более удачным образом.

Вывод

Буду откровенен – работает оно так себе. Все графики и прочие результаты я продемонстрировал, можно видеть, что результат работы оставляет желать лучшего. Во всяком случае, для поиска точного решения на данном этапе этот алгоритм использовать нельзя. С другой стороны, он все-таки находит более-менее вменяемые варианты, в откровенную лажу не скатывается. Если задача не требует точного решения, а сойдет и близкое к оптимальному, то этот алгоритм вполне себе работоспособен. Как вариант, можно использовать этот алгоритм для начального приближения к оптимальному значению, а затем запускать что-нибудь типа перебора вариантов.
Если у вас будут какие-то идеи улучшения работы алгоритма, буду очень рад их выслушать.
Также готов выслушать любую конструктивную критику по самой статье (оформление, стиль изложения, etc.).

Для желающих покопаться в кодах самим выложил на mail cloud архив с проектом (написано в VS2012, с более ранними версиями VS проект будет несовместим; для запуска программы потребуется .net framework 4.5, .exe-файл находится с папке \Salesman\Salesman\bin\Debug)