CC BY
54
УДК 519.6
ИСПОЛЬЗОВАНИЕ ПСЕВДОСЛУЧАЙНЫХ ВЕЛИЧИН В МЕТОДЕ ГРАВИТАЦИОННОЙ КИНЕМАТИКИ ПОИСКА ГЛОБАЛЬНОГО ЭКСТРЕМУМА
А.В. ПАНТЕЛЕЕВ, Е.А. АЛЁШИНА
Рассмотрены модификации метода гравитационной кинематики с использованием псевдослучайных величин, позволяющие находить глобальный экстремум функций со сложной структурой поверхностей уровня. Разработано программное обеспечение метода для решения задач условной оптимизации с применением рассмотренных модификаций. На тестовых примерах исследуется эффективность модификаций метода, а также возможность их применения к решению практических задач.
Ключевые слова: псевдослучайные величины, метод гравитационной кинематики, глобальный экстремум.
Введение
Метод гравитационной кинематики, предложенный R.A.Formato в[1], относится к группе ме-таэвристических методов оптимизации. Использование таких методов становится актуальным для решения разнообразных задач оптимального управления летальными аппаратами ввиду большой размерности решаемых задач, а также сложной структуры математических моделей.
Рассматриваются модификации метода гравитационной кинематики, использующие псевдослучайные величины. Основное свойство этих модификаций состоит в том, что они полностью сохраняют детерминированность самого метода, что обеспечивает воспроизводимость экспериментов.
Авторами статьи разработано алгоритмическое и программное обеспечение для решения задач условной оптимизации с применением модификаций метода. Решение тестовых примеров показывает, что использование модификаций с псевдослучайными величинами улучшает функционирование метода гравитационной кинематики и позволяет ему конкурировать с методом частиц в стае или непрерывной модификацией метода муравьиных колоний, рассмотренными авторами ранее [2, 3].
1. Постановка задачи
Дана непрерывная функция f (х) = f (х1, х2,..., xn), определенная на множестве допустимых решений D с Rn, где х = (х1, х2,..., xn)T, D = {х | xt е [at,bt], i = 1,2,...,n} .
Требуется найти условный глобальный минимум функции f (х) на множестве D, т.е. такую точку х* е D, что f (х*) = min f (х).
xеD
2. Стратегия поиска решения
Идея метода гравитационной кинематики подробно описана в [1,4]. Рассматривается эволюция популяции, содержащей m точек (частиц): х1,..., хт . Метод использует аналогию с законом всемирного тяготения. Считается, что каждая частица популяции движется в множестве допустимых решений D под действием гравитационных сил. Аналогом массы является приращение целевой функции.
Частица с номером ] из положения х]’к переходит в новое положение х]’к+1 по траектории, определяемой начальным положением, скоростью и ускорением, порождаемым силой взаимодействия с остальными частицами.
Ускорение у -й частицы в силу взаимодействия с частицами, имеющими номера 1,2,...,] -1,] +1,...,т , находится по формуле
’ ,к / ,к
а« = С-21 [/(х“)-/(х‘-‘)]■[/(х',к)-/(х‘-к)]" ■^-Х-ь; (1)
’*- - \хик - х],к
1 [. 1 (х’")-/(х',к)'
где С - аналог гравитационной постоянной У; множитель [/(х’,к)-/(х1,к) (приращение
значения целевой функции) - аналог массы; множитель (х1,к - х1,к) - аналог вектора Г в законе
I ’ к 1 к \Р 3 П
всемирного тяготения; выражение х’ -х^ - аналог величины г ; а, р - параметры, гаран-
( [1, * > о ^
тирующие гибкость поиска; 1 - единичная ступенчатая функция 1(*) = < ; выражение
I I0, * <0;
введено с целью не допустить появления отрицательных масс.
Скорость частицы для упрощения расчетов можно положить равной нулю, считая частицы в каждый момент времени неподвижными.
Тогда будущее положение частицы находится по следующей формуле
хМ+1 = х,к + - аьк М2. (2)
2
После нахождения новых пробных точек производится проверка их принадлежности области допустимых решений (ОДР). Если точка не принадлежит ОДР, необходимо принять меры по
пересчету положения. Далее либо среди всех точек {х1к,..., хт,к, х1,к+1,..., хт,к+1} выбирается т наилучших по величине функции, либо в качестве новой популяции используется множество
{1 к+1 т к+1 "1 ^
х , . , х } новых положений точек.
Процедура перехода от одной популяции к другой заканчивается по достижении заданного их числа К.
3. Модификации с использованием псевдослучайных величин
При сравнении метода гравитационной кинематики с ранее рассмотренными авторами методами муравьиных колоний и частиц в стае [2, 3] видно, что случайный характер других метаэври-стических методов все же дает им преимущество перед методом гравитационной кинематики [4]. В [5] была предпринята попытка внедрить псевдослучайные элементы. Такой способ улучшает характеристики метода, оставляя его детерминированным. Под «псевдослучайностью» понимается использование детерминированной, но меняющейся со временем, величины независимо от самого процесса поиска, от структуры ОДР и целевой функции. Псевдослучайный подход обеспечивает полную воспроизводимость работы алгоритма, т.е. сохраняет отличительную черту исходного метода гравитационной кинематики. При этом модифицированный алгоритм должен выдавать лучшие результаты по сравнению с обычным за счет «встряхивания» координат частиц.
В [4] рассматриваются три способа внедрения псевдослучайных величин: использование псевдослучайного начального распределения частиц в ОДР; применение псевдослучайного множителя при принудительном возвращении вылетевшей за пределы ОДР частицы; сжатие пространства поиска вокруг лучшей на рассматриваемом этапе позиции.
Рассмотрим каждый из способов подробнее.
Модификация 1. Использование псевдослучайного начального распределения частиц, которые располагаются равномерно вдоль пересекающихся прямых, параллельных координатным осям. Схематическое изображение распределения представлено на рис. 1. Точка пересечения прямых находится на диагонали ОДР (п -мерного параллелепипеда), ее положение задается через радиус-вектор ё, компоненты которого находятся по следующей формуле
di = а + у(Ь - а ), г = 1, п . (3)
Параметр /е [0,1] является псевдослучайным. Его значение задается при каждом запуске программы неслучайным образом, описанным далее.
Рис. 2. Сжатие границ ОДР вокруг лучшей на данной итерации пробной точки
Рис. 1. Псевдослучайное начальное распределение частиц
Модификация 2. При выходе частицы за границы ОДР ее новое положение рассчитывает-
ся по следующей формуле: если х/,к < а, то полагают
Xм
,м
і ' аг + РгеР (ХҐ - аі ) , і = 1П , (4)
если х1і’к > Ьг, то
ХҐ = Ь,+ рг„ (Ь - х‘*), г = Щ, (5)
где Егер - псевдослучайная величина, ее значение меняется по закону, не связанному ни
со структурой ОДР, ни с целевой функцией.
Модификация 3. Сжатие ОДР вокруг лучшего на данный момент решения через определенное число итераций Т. Схема процесса представлена на рис. 2. Новые границы пространства поиска задаются следующим образом
Ь - хи —
Ь Х • - (6)
а = а +-
хи - а,.
Ь = Ь
і = 1, п,
2 - ‘ ‘ 2
где хи - позиция частицы, соответствующая наименьшему значению функции, полученному за количество итераций, кратное Т. Вектор хи не является случайным, так как ему не ставится в соответствие какой-либо закон распределения. После процесса сжатия часть точек может оказаться вне ОДР. Для этих точек можно сгенерировать новые позиции, которые будут удовлетворять новым ограничениям, или же не менять их положения, приняв во внимание, что в новую ОДР искомая точка минимума может и не попасть.
4. Алгоритм с учетом модификаций
При изложении алгоритма используется покоординатная форма записи.
Шаг 1. Задать К - максимальное число итераций; т - число пробных точек в текущем поколении; параметры метода G , а , Ь, At. Задать, если требуется, параметры модификаций /, Т. Шаг 2. Положить к = 0 .
псевдослучайным образом:
а) если не используется псевдослучайное начальное распределение точек, расположить их в виде сетки на ОДР;
б) если выбрана Модификация 1, задать значение параметра /, найти радиус-вектор точки пересечения прямых, на которых будут располагаться пробные точки, по формуле (3), а затем по вектору ё сгенерировать начальные позиции.
Шаг 2.2. Подсчитать значения целевой функции / (х10),..., / (хт 0).
Шаг 3. Генерирование новых пробных точек.
Шаг 3.1. Положить ] = 1.
Шаг 3.2. Пересчитать значения компонент ускорения для ] -й точки популяции в соответствии с формулой (1).
Шаг 3.3. Найти новое положение ] -й точки популяции, используя формулу (2).
а) если Модификация 2 не используется, то величина Dt в (2) уменьшается до тех пор, пока не выполнится условие x.’k+1 е [,bt];
б) если выбрана Модификация 2, то новая позиция находится по формулам (4), (5).
Шаг 3.4. Если j = m, процесс завершить и перейти к Шагу 4. Если j < m, то положить j = j +1 и перейти к Шагу 3.2.
Шаг 4. Формирование новой популяции.
Шаг 4.1. Сформировать множество {xu,...,xm,k,x1,k+1,...,xm,k+1} ; подсчитать значения функции в пробных точках f (xu+1),..., f (xm,k+1).
Шаг 4.2. Отобрать m наилучших решений, а остальные отбросить.
Шаг 5. Если k < K, то продолжать поиск решения:
а) если выбрана Модификация 3, и k mod T = 0, задать новые границы ОДР по формулам (6);
б) положить k = k +1 и перейти к Шагу 3.
Если k = K, процесс закончить, в качестве приближенного решения выбрать наилучшую точку из последней популяции.
5. Программное обеспечение
На основе разработанного алгоритма сформировано программное обеспечение на языке C# в среде Microsoft Visual Studio 2005. Пользовательский интерфейс аналогичен интерфейсу программы, разработанной для изучения методов муравьиных колоний и частиц в стае и рассмотренной в [2, 3]. С помощью программного обеспечения пользователь может выполнять следующие действия:
- вводить параметры постановки задачи;
- задавать параметры метода;
- наблюдать процесс поиска решения в виде анимации;
- использовать модификации алгоритма;
- анализировать полученный результат;
- сохранять результат в памяти компьютера для последующего анализа.
Шаг 2.1. Генерировать m точек в множестве
детерминированным или
Если xj,k+1 £ [at, b ], то
На рис. 3 показана структура интерфейса окна выбора модификаций метода гравитационной кинематики.
Рис. 3. Окно выбора модификаций метода гравитационной кинематики 6. Решение тестовых примеров
Для исследования модификаций метода гравитационной кинематики рассмотрим задачи минимизации двух овражных функций, Розенброка и Three Humps Camel Back и двух функций Акли и Швефеля, с множеством локальных минимумов (табл. 1, рис. 4).
Таблица 1
Название Формула ОДР * x f *
Функция Розенброка x2) = 100 (x2 - X2 )2 +(1 - )2 -5 < x1, x2 < 5 (1,1) 0
Функция Three Humps Camel Back x6 f2 (X1 , X2) = 2X12 - 1 05X14 +4- + X1X2 + X22 6 -5 < x1, x2 < 5 (0,0) 0
Функция Акли f3(x1, x2) = 20 1 - exp (-0,2^0,5(x2 + y2)) -- exp (0,5 [cos(2^x) + cos(2^y)j) + e -35 < x1, x2 < 35 (0,0) 0
Функция Швефеля f4( ^.^ x„ ) =-Ydxi ) i=1 -500 < x < 500 * x = = 420,96 -418.98 •n
Рис. 4. Графики функций: а - Розенброка; б - Three Humps Camel Back; в - Акли; г - Швефеля в окрестности точки глобального минимума (а, б) и на всей ОДР (в, г)
Сравним различные модификации с псевдослучайными величинами при следующих значениях параметров: m = 180; G = 2; a = 2; ¡3 = 2; K = 500; At = 1. При использовании модификации 1 значение параметра g варьировалось от 0 до 1 с шагом 0,1; модификация 2 использовалась в соответствии с рекомендациями в [5], в модификации 3 параметр T принимал значения 10, 20, 30, 40. Далее из всех результатов выбирались минимальные значения, которые были сгруппированы в табл. 2.
Таблица 2
Модификации Минимальные значения из всех полученных
1 2 3 А f2 f3 a = 2
- - - 0,006975277 0,002245527 0,068502327 -837,906177
+ - - 0,009665762 3,94264E-05 2,28144E-10 -837,965760
- + - 0,006975277 0,003120362 0,068502327 -837,965774
- - + 1,06966E-05 2,22784E-09 0,000167916 -837,965775
+ + - 0,012623669 4,94684E-07 2,28144E-10 -837,965581
+ - + 0,009665762 8,14015E-07 4,76239E-10 -837,965500
- + + 0,000196469 6,11342E-10 0,000220068 -837,965775
+ + + 0,000151963 4,9686E-09 4,76239E-10 -837,965773
По значениям в табл. 2 заключаем, что наиболее эффективна модификация 3, а также сочетание модификаций 2 и 3, или использование всех трех модификаций одновременно. Для функции Акли эффективной также оказалась модификация 1. Не дало улучшения результатов только использование одной модификации 2 (ее лучше применять вместе с модификациями 1 или 3).
В табл. 3 приведены результаты решения задач минимизации функций Акли (/3* = 0) и
Швефеля при п = 2 (/4 » -837,97 ), при различных значениях параметра у от 0,1 до 0,9 с шагом 0,2 при т = 120, т = 180 . Решение задачи производилось при следующих значениях параметров: О = 2; а = 2; /3 = 2; К = 500 ; At = 1; Т = 20 для /3 и Т = 40 для /4 .
Таблица 3
Значения параметров Значение функции /3 Значение функции /4, п = 2
у т О а 3 К At
0,1 120 2 2 2 500 1 19,58495145 -719,5273999
0,3 120 2 2 2 500 1 0,00136526 -719,5272723
0,5 120 2 2 2 500 1 0,02689170 -697,7688620
0,7 120 2 2 2 500 1 1,38092Е-10 -809,2096720
0,9 120 2 2 2 500 1 19,56069681 -837,9655233
0,1 180 2 2 2 500 1 18,00445516 -719,5274283
0,3 180 2 2 2 500 1 0,17165491 -696,3705563
0,5 180 2 2 2 500 1 4,76239E-10 -837,5906997
0,7 180 2 2 2 500 1 1,93505E-09 -837,9657730
0,9 180 2 2 2 500 1 19,25041664 -837,9652203
Наконец, рассмотрим задачу большой размерности (п = 30). Будем минимизировать функцию Швефеля (/4 »12569,49). Результаты решения задачи при различных значениях псевдослучайных параметров приведены в табл. 4.
Таблица 4
Значения параметров Значение функции /4, п = 30 Отклонение от точного решения
у т О а 3 К At Т = 40
0,5 120 2 2 2 500 1 20 -3671,31 8898,18
0,7 120 5 5 2 500 1 20 -9575,37 2994,12
0,7 120 2 2 2 500 1 20 -5980,52 6588,97
0,9 120 2 2 2 500 1 40 -11222,68 1346,81
0,9 120 2 2 2 500 0,5 40 -12216,16 353,33
0,95 120 2 2 2 500 1 40 -12128,85 440,64
0,95 120 2 2 2 500 0,5 40 -12289,39 280,10
1,0 120 2 2 2 500 1 40 -10902,42 1667,07
0,5 180 2 2 2 500 1 40 -5699,36 6870,13
0,7 180 2 2 2 500 1 40 -11981,11 588,38
0,9 180 2 2 2 500 1 40 -10889,90 1679,59
0,95 180 2 2 2 500 1 40 -11841,63 727,86
Учитывая сложную структуру функции и большую размерность задачи, полученные результаты можно считать приемлемыми. Улучшить качество решений можно путем дальнейшей настройки параметров, а также комбинированием метода гравитационной кинематики с другими метаэвристическими методами оптимизации.
Заключение
Рассмотрены различные способы улучшения метода гравитационной кинематики за счет использования «псевдослучайных» параметров. Для исследования модификаций метода с «псевдослучайными величинами» разработано программное обеспечение. Показано, что для функций сложной структуры модифицированный метод приводит к нахождению в целом более точных решений, чем метод гравитационной кинематики без модификаций. Модифицированный метод гравитационной кинематики способен конкурировать и с другими метаэвристиче-скими методами оптимизации. При этом использование именно псевдослучайных, а не случайных параметров обеспечивает воспроизводимость процесса решения задачи.
ЛИТЕРАТУРА
1. Formato, R.A. Central Force Optimization: а new metaheuristic with applications in applied electromagnetics, Progress In Electromagnetics Research, PIER 77, 425-491, 2007.
2. Пантелеев А.В., Алёшина Е.А. Метод частиц в стае // Метаэвристические алгоритмы поиска глобального экстремума. - М.: Изд-во МАИ-ПРИНТ, 2009. - Гл. 4. - С. 85-106.
3. Алёшина Е.А. Применение метода муравьиных колоний к задаче поиска условного глобального минимума функции многих переменных // Теоретические вопросы вычислительной техники и программного обеспечения: межвуз. сб. научн. тр. - М.: МИРЭА, 2011. - С. 63-68.
4. Пантелеев А.В., Алёшина Е.А. Разработка алгоритмического и программного обеспечения метода гравитационной кинематики // Научный Вестник МГТУ ГА. - 2011. - № 165. - С. 95-102.
5. Formato R.A. Pseudorandomness in Central Force Optimization // Статья на сайте www.arXiv.org.
USE OF PSEUDORANDOM VARIABLES IN CENTRAL FORCE OPTIMIZATION
Panteleyev A.V., Alyoshina E.A.
We consider a process of searching the conditional global extremum of multivariable function. The solution is searched by using modified Central Force Optimization (CFO). The modifications include pseudorandom variables, this is significant because the whole process is then reproducible opposite to other metaheuristics. The developed software proves that that approach is efficient; this is shown on some test examples including multidimensional functions with many local minima and functions with sophisticated structure. This allows to apply modified CFO to real problems in future.
Key words: pseudo-random values, method of gravitational kinematics, global extremum.
Сведения об авторах
Пантелеев Андрей Владимирович, 1955 г.р., окончил МВТУ им. Н.Э. Баумана (1978), доктор физико-математических наук, профессор, заведующий кафедрой математической кибернетики МАИ, автор более 150 научных работ, область научных интересов - методы синтеза оптимальных нелинейных систем управления, методы оптимизации.
Алёшина Екатерина Александровна, окончила МАИ (2010), аспирантка МАИ, автор 10 научных работ, область научных интересов - методы оптимизации, численные методы.
